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Abstract 


A computer program for evaluating the electromagnetic field pattern of a 
known arbitrary incident field scattered from a perfectly conducting reflector of 
arbitrary shape is presented. 

It is shown that the commonly used assumption of far-field behavior for the 
incident field leads to poor results in some cases; an exact method is developed in 
which a spherical-wave expansion is used to represent the incident field, and it is 
shown that this method yields the correct result in these cases. 

The scattering surface is represented by a Fourier type of expansion in one 
variable, where the coefficients are specified as tabular functions of the remaining 
variable. A tilted figure of revolution is a case of practical importance, and it is 
shown that, for small tilt angles, three Fourier components are usually sufficient 
to represent the surface. 

A new nonlinear integration technique is used that is four to eight times faster 
than Simpson’s rule integration, under specified conditions. Other techniques to 
maximize program efficiency are described, and a comparison is given that indi- 
cates the program to be 11 to 19 times faster than a similar existing program. 

The subreflector and main reflector patterns of an 85-ft-diam Cassegrainian 
antenna are computed, using both the far-field approximation and a spherical- 
wave expansion to represent the incident fields. The results are compared, and it 
is shown that the spherical-wave expansion of the subreflector pattern agrees in 
remarkable detail with an experimentally measured pattern. 

Exemplifying the general capabilities of the program, the patterns of a para- 
boloid with a displaced feed are computed for three cases, and the results are 
compared with approximate theory. 
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Calculation of Scattered Patterns 
From Asymmetrical Reflectors 


I. Introduction 

The physical-optics technique has been used with great 
success to compute the scattering from perfectly con- 
ducting surfaces, such as the subreflector of a Casse- 
grainian antenna (Ref. 1). In fact, the technique has 
yielded such accurate results that computed patterns 
have begun to replace a good deal of experimental data 
in antenna development and analysis (Refs. 2 and 3). 

A major limitation of the physical-optics technique is 
that it requires the evaluation of a double integral with 
a rapidly oscillating integrand that typically results in 
the need for large blocks of stored data and increased 
computer time, In some cases, with special symmetry 
conditions, one integration may be performed analyti- 
cally (Ref. 4). The program described in this report was 
developed so that the broad class of scattering problems, 
which do not necessarily have any symmetries, could be 
resolved in the most efficient manner possible. The great 
majority of computer time and storage requirements may 
be directly attributed to the numerical evaluation of the 
integral; therefore, the major contribution towards effi- 
ciency in this program is the application of a new 
integration technique (Ref. 5). 


The case of scattering from an infinite plane reflector 
is one of the few electromagnetic scattering problems 
that has a simple known solution; the resulting field is, of 
course, identically the reflection of the incident field. 
An earlier version of the computer program described 
herein implicitly assumed far-field behavior for the 
incident-field pattern. When the program was checked, 
with the plane reflector problem as a test case, it was 
found that the computed scattered pattern was not the 
same as the incident pattern and that the result was 
dependent upon the wavelength of the incident field. 
(Data for this test case are presented in detail in 
Section III.) 

Of the several possible explanations for this anomaly, 
it was hypothesized that the far-field behavior of the 
incident fields was the cause. A technique based upon 
spherical-wave theory (which yields an exact repre- 
sentation of the incident fields) was developed, and the 
application of this technique eliminated the error. 

The spherical- wave representation of the incident fields 
is in principle completely general. Also, a general repre- 
sentation of the scattering surface was used in the 
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development of the program. Practical considerations 
required that some of this generality be sacrificed in 
the actual programming, but the theory could easily be 
applied to extend the results if required. 

The only important assumption made in the develop- 
ment of this program was the physical-optics approxi- 
mation in which the currents on the reflecting surface 
were assumed to be determined by the incident magnetic 
field. The validity of this procedure is discussed in 
Section II; however, the agreement between experimental 
and computed results presented in Section VI and else- 
where (see Refs. 1-3) is, in fact, the strongest demon- 
stration that the physical-optics approximation is a useful 
and valid engineering method. 

II. The Physical-Optics Technique 

If the currents K induced on a surface S by incident 
fields E i and H * can be determined, then the scattered 
fields E s and H s caused by these currents are given by 

Es(P) = ~ldb / [(K. • V) V + k 2 K] dS 

(la) 

and 

Hs(p) = i!r/ S (KX v)£ r- ds ( lb ) 

The coordinate system is shown in Fig. 1. 

Since these equations (Ref. 6) are exact, the scattering 
problem consists of (1) determining the currents K and 
(2) evaluating the integrals. 

One method for obtaining the currents is to use the 
following boundary condition on S: 

n X (Ej + E a ) = 0 (2) 

This method leads to the following integral equation for 
the currents &; 

— n X Ej = — n X -4 — [ [(K* V) V + fc 2 K] dS 

47ra>€ J s T 

(3) 

This formulation has been known for some time 
(Ref. 6), but the speed and size of the present generation 
of computers have made the technique much more 


attractive. This integral equation method has been used 
recently with great success in solving two-dimensional 
scattering problems and problems involving a surface of 
revolution (Ref. 7). The major disadvantage of this 
approach is that, in these restricted cases, the surface 
must not be wider than approximately 20 wavelengths 
(Refs. 8 and 9). In the case of an arbitrary surface, the 
limitation is more severe. 

An approximate alternate method for obtaining the 
currents is to assume that, on areas directly illuminated 
by the source, the currents are the same as they would 
be if the incident fields were reflected optically: 

K — 2n X (4) 

On shadowed regions, the currents are assumed to be 
zero. These assumptions are known as the physical-optics 
approximation. 

For smooth reflectors, this approximation is valid in the 
limit of zero wavelength (in Ref. 10, Cullen gives proof 
of this for the case of a convex body illuminated by a 
plane wave) and, in contrast to the integral-equation 
method, the approximation improves as surface 
dimensions become large compared to a wavelength. 

However, not only does this approximation yield 
results when it is not practical to apply the integral- 
equation technique, but it yields the results almost free 
in terms of computer time. Therefore, the generally 
time-consuming integral-equation technique should be 
used only when the physical-optics approximation is 
inadequate. In the typical scattering problem, the cur- 
rents are of little direct interest; they are only an inter- 
mediate step toward obtaining the scattered fields. 
Therefore, the physical-optics approximation should not 
be judged by whether it yields the correct currents, but 
by whether the final computed scattered fields are correct. 

In some cases Eq. (4) is a poor approximation of the 
true currents, but it results in an excellent approximation 
of the true scattered fields. The reason for this may be 
understood by the following: let the true currents be K, 
and define AK by 

aK = K - 2nXHj (5) 

Typically, AK exhibits oscillatory behavior over the 
surface, and its net integrated contribution to the scattered 
fields is small. Watson (Ref. 11) has shown this analyti- 
cally for the particular case of the fields in the focal 
region of a paraboloid that is large with respect to a 
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wavelength. Presumably, the primary effect of AK is 
to produce local fields that represent stored rather than 
propagating energy. 

The similarity between the currents computed by the 
integral equation technique and the physical-optics 
approximation has also been pointed out by Andreasen. 
The oscillatory nature of the difference may be seen in 
his results that are for scattering surfaces only a few 
wavelengths wide (Ref. 9). Surface currents for small 
spheres — 0.18 to 3.2 wavelengths radius — obtained from 
a classical boundary value solution are given by King 
and Wu in Ref. 12, and even here, Eq. (4) is a reasonable 
approximation. 

Therefore, Eq. (4) is not only a necessary step toward 
reaching a result in problems involving a large reflector, 
but also a sound engineering approximation for surfaces 
larger than a few wavelengths in size. 

If it is assumed that n and Hi are known, the integrals 
of Eq. (1) can be evaluated to obtain the scattered fields 
at any point. However, the problem is greatly reduced if 
it is assumed that E s and H s are to be evaluated only at 


very large distances from S. As R-» oo, E s and H s must 
satisfy the far-field relations; 

H = (^-) V i« XE (6a) 

and 

z ?— jJcR 

E(R,®,4>) - [Ee(0,*) ie + £*(©,$) h] (6b) 

Therefore, only E© and E$ need be evaluated. Also 
in this case, Eq. (1) may be written so that the only term 
in the integrand that is dependent on the location of the 
field point P is the argument of the exponential term. 
The remainder of the integrand may be initially com-' 
puted over S and stored to be used repeatedly in the 
evaluation of the integral for every field point P. This 
method represents a substantial reduction in computer 
time. 

The assumption that the scattered fields are to be 
evaluated at R = oo does not represent a loss of gen- 
erality. If near-field results are required, they may be 
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obtained from the far-field data using the spherical-wave 
expansions described later. 

The far-field version of Eq. (1), with the currents that 
are obtained from Eq. (14), is derived in Silver (Ref. 6): 

p-jkR 

E.(P) - Ieie + hU) (7) 


where 


I = — J (n X Hi) exp (jkp * i I{ ) dS 

The final result, the total field E r , is obtained by adding 
the incident fields to the scattered fields as shown in 
Eq. (8): 


E T — E s + E* (8) 

The problem has now been reduced to the following 
steps: (1) specify the surface data p and n, (2) specify 
the incident field data H f , and (3) evaluate the integral. 


III. Specification ©f the Scattering Surface 

In the derivation of the scattering integral equations 
(see Eq. 1), the surface S is assumed to completely en- 
close some region of space. Examples of surfaces are a 
sphere (which separates space into an infinite region 
and a finite region) and an infinite plane (which sep- 
arates space into two infinite regions). There are obvious 
numerical difficulties with infinite surfaces; the finite 
closed surface also presents some potential problems 
which are discussed in the following paragraphs. 

As discussed in Section IV, the origin of the coordinate 
system is determined by the location of the source, and 
the natural system in which to describe the incident 
fields is spherical coordinates. To define a surface, it is 
necessary and sufficient to specify one variable as a 
function of the other two. The most natural and con- 
venient choice for this problem is a function p(0,<f>). In 
the case of a sphere that does not enclose the origin, a 
line 0 = const <f> = const, in general, intersects the surface 
twice so it cannot be described by a single function of 
this form. However, with the assumed current approx- 
imation, the back portion of the sphere has zero currents 
and can be ignored. The illuminated portion is precisely 
the portion that can be defined by a single function. 
There is still one potential difficulty. At the point at 


which p is just tangent to the surface, the partial deriva- 
tive dp/dO is infinite and, as will be seen shortly, is re- 
quired. Therefore, it is necessary to discard a further 
portion of the surface. In the case of an infinite surface, 
all but some finite portion is ignored. In both cases the 
remainder may be identified as a truncated surface. 
Therefore, an ‘arbitrary reflecting surface” is defined as 
a surface that can be represented by a function p(0,<l>) 
whose partial derivatives exist at every point on the 
surface. This definition is basically a consequence of 
the physical-optics approximation used to obtain the 
currents, and is sufficiently general to include virtually 
all antenna reflectors. (Discontinuous surfaces may be 
represented by two or more segments satisfying these 
requirements, and the results may be superimposed.) 


Because this surface may be thought of as a portion 
of a closed surface, the derivation of Eq. (1) is still valid, 
except for one point. A truncated surface has an edge, 
and the currents obtained from Eq. (14) will in general 
be discontinuous there. To ensure that the results are 
consistent with Maxwell’s equations and the radiation 
condition, Silver (see Ref. 6) introduced a line charge 
on this edge in the derivation of Eq. (7). Sancer (Ref. 13) 
has shown that the surface integrals of Eq. (1) intrinsi- 
cally contain the contour integral representing the con- 
tribution of this charge, and that the term introduced by 
Silver arises automatically as a mathematical consequence 
of truncating the surface. 


The surface data required to evaluate the integral of 
Eq. (7) are p and n. Using a result from differential 
geometry (Ref. 14), Dickinson (Ref. 15) pointed out the 
useful relationship : 


ndS = X — 

0<p ou 


where 


P = pi P 

dp dp . . . 

— =- lp + P sme H 


dp _ 0p . 

lo~ d0 lp 


pio 


(9) 


( 10 ) 


Therefore, specification of p and its partial derivatives 
is a complete description of the surface. Since these 
quantities must be specified on a set of points on S, which 
will be called the integration grid and can easily involve 
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several thousand points, the storage and transmission of 
a description of a surface in this form can be a problem. 
For this reason, it is useful to achieve some sort of data 
compression. One trivial case is that in which the surface 
is completely described by an analytic expression; for 
example, a paraboloid whose focus is at the origin, is 
completely described by a single number f as shown in 

Eq* (11): 


P=f sec 2 — (11) 

A compromise between the complete generality of two- 
dimensional tabular input and the severe restriction of a 
simple analytic equation may be developed from the 
general form : 

M 

p(0,<f>) = a m (0) cos m<f> + b m (0) sin m<f> (12) 

m~ o 

Then 

M 

dp 0a m . , , • , /TO \ 

w = cos + ~W sm (13a) 

•m=o 


and 


M 

Z/ 


— mu m (0) sin m<f) + mb m (&) cos m<t> 


(13b) 


Suppose that for a given surface it is necessary to 
include terms of order up to M to determine p with 
sufficient accuracy; then by the sampling theorem it 
would be necessary to specify p at 2 M values of <j> to 
accurately represent details of the surface variations using 
tabular data. Therefore, this Fourier type of represen- 
tation is in principle always at least as efficient as an 
ordinary tabular representation. 


If M is allowed to be arbitrarily large, Eq. (12) is in 
fact completely general; however, the reason for choosing 
this form is that the first few terms are sufficient to 
accurately represent a slightly tilted surface of revolution, 
which is a case of substantial practical importance. 
For example, Fig. 2 shows the peak error that occurs 
when a tilted plane reflector is represented by three terms 
of the expansion, a o (0 ), ajjf), and a 2 (0). (The axis of 
rotation is taken to be the y-axis, so all of the b m (Q) terms 
are zero.) This error is in general a function of the curva- 
ture of the reflector, but the data shown in the figure are 


o 

z 

< 



0 20 40 60 

MAXIMUM ANGLE 8 SUBTENDED BY REFLECTOR, deg 
max 


Fig. 2. Peak error with three Fourier components used to 
represent a tilted reflector 


representative of the accuracy that may be obtained 
with this method. A reasonable a posteriori estimate of 
accuracy may be obtained by extrapolating the value of 
the first neglected term in the expansion. For example, 
for a tilt angle of 1 deg, the coefficients at 0 ma x = 20 deg 
are a 0 = 1.064199, a, - 0.006761, and a 2 = 0.000021. The 
values decrease with a nearly constant ratio, and a 
reasonable estimate for the neglected third-order term is 
about 10~ 5 to 10~A which is approximately the peak error. 
The relatively poor accuracy obtained for a tilt angle of 
10 deg at # nmx = 40 deg is also predictable from the co- 
efficients (a 0 = 1.320015, a ± = 0.197465, and a 2 = 0.014608) 
and would lead to an expected relative error of approxi- 
mately 10~ 3 . 

The derivatives of the a m (6) terms must also be written 
in as tabular functions of 8; therefore, six input functions 
are required in this example where M= 3 and all b m (0) 
terms are identically zero. The number of cf> values at 
which p must be known is determined by the require- 
ments of the numerical integration technique and is 
usually quite large (from fifty to several hundred values 
is typical). Therefore, the specification of only six func- 
tions is a substantial reduction in the quantity of data 
required. A computer program to obtain these functions 
for the case of an arbitrary tilt and translation of a surface 
of revolution is described in Appendix C. 


IV. Specification of the Incident Fields 

The most common method (Refs. 1, 6, and 16) for 
obtaining the magnetic field on the surface S is to assume 
that the incident fields satisfy the far-field relations: 


H,- = 



(14a) 
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and 


g-i&P 

Ei ( p , 0, <f>) = [Eg ( 0, <£) i e + E^ (0, <£) 4] — — (14b) 

Not only is this a relatively simple form, but Ee (0, 0) 
and E^ (0, </>) correspond to the quantities conventionally 
measured on an antenna pattern range, so experimental 
data may be easily used as input for the incident fields. 
Also, as discussed in Section II, theoretical patterns are 
easier to compute in this form. 

Equation (14) represents a spherical wave diverging 
from the origin; therefore, it is not valid unless the origin 
is located at the source phase center. The phase pattern 
about this point should be constant, but in practice the 
relative phase may vary by 90 deg or more over the 
region of interest. A nonconstant phase pattern may be 
taken into consideration if E e and E^ are allowed to be 
complex valued, but the radial components of the field 
implied by the nonconstant phase pattern are usually 
assumed negligible. 

The use of this assumption provides good results for 
many practical problems. Serious errors can arise, how- 
ever, when the surface is “close” to the source. This fact 
has been recognized by Zucker and Ierley, who directly 
computed the near-field radiation of a conical horn in eval- 
uating the scattered pattern of a near-field Cassegrainian 
subreflector (Ref. 3). An alternative approach is to take 
experimental data at range values corresponding to the 
location of the surface; this was done by Hogg and 
Semplak (Ref. 17). 

The disadvantage of the first method is that horn 
patterns computed in this way frequently do not agree 
with experimental data, since it is necessary that strong 
assumptions be made about the fields in the aperture of 
the horn. The disadvantage of the second method is that 
each different subreflector configuration in principle 
requires a new set of experimental data. These problems 
are obviated by the use of the spherical-wave expansions 
that are discussed later in this section. 

A graphic demonstration of the error that may be intro- 
duced if far-field behavior is assumed for the incident 
fields is provided by the case of scattering from a large 
plane reflector. In the case of an infinite plane, the result 
must be the reflection of the incident pattern, inde- 
pendent of frequency. Also, in the case of an infinite 
plane, Eq. (4) is valid, independent of frequency; there- 
fore, this is an excellent direct test of Eq. (14). 


The cases shown in Fig. 3 (and one additional case at 
z — 160 A, which is not illustrated) were evaluated using 
the far-field form for the incident fields given by Eq. (14). 
The fields were derived from the experimental pattern 
of a conical horn with an aperture diameter of 4.671 A, 
illustrated at the origin of Fig, 3. The actual numerical 
integration was truncated at an angle at which contri- 
butions became negligible. The disks shown in Fig. 3 
represent the region over which the integral was actually 
evaluated. It should be noted that horn-reflector inter- 
actions are not considered, so this is not intended to be 
a perfect duplication of a real experimental situation, but 
rather a numerical test of the far-field assumption given 
by Eq. (14). 

The computed scattered fields based on the far-field 
approximation of the incident field for the six cases are 
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shown in Fig. 4 along with the incident-field pattern. 
Only the ^-component is illustrated, but the behavior of 
the ©-component is similar. All of these patterns should 
be identical, but the discrepancies are large. In addition 
to the severe distortion of the reflected main lobe when 
z — 5 and 10 A, a spurious back lobe appears in all cases 
— that is, the plane that was assumed to be a perfect 
reflector behaves as if it were partially transparent. 

In the special case of a plane reflector, the scattered 
fields E s are identical in the reflected and backward 
directions. When the incident field is added (Eq. 8), it 
should cancel this backward radiation. Therefore, this 


back lobe is identically the error between the reflected 
pattern and the incident pattern. (The error includes 
both amplitude and phase differences since it is the 
difference between complex-valued quantities.) 

The power in the computed scattered patterns has been 
calculated as a percentage of the power in the incident 
fields and is shown in Fig. 5. As indicated in Fig. 5, the 
main (reflected) lobe alone accounts for 100 percent of 
the incident power in all cases, so the total power in the 
computed pattern exceeds the incident power. This excess 
is, of course, physically impossible. The spurious power 
in the back lobe is strongly dependent on the distance 



0 36 72 108 144 180 

POLAR ANGIE 8, deg 

Fig. 4. Computed amplitude patterns for plane reflector cases; far-field assumption for incident fields 
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from the source to the surface, which is plotted in units of 
D 2 /A where D is the diameter of the source. At £=2D 2 /A, 
the traditional far-field boundary line, this power has 
dropped to less than 2 percent. 

The behavior of the phase pattern of E$ is shown in 
Fig. 6. These phase data have been transformed to a new 
coordinate system; the origin was moved from the point 
at the source phase center to the image of this point 
behind the plane reflector where the phase center of the 
total scattered fields should be located. In fact, this is 
very nearly the phase center for the first three cases; the 
patterns of the last cases are those of a source whose 
phase center is closer to the plane surface. 

Therefore, instead of the total scattered fields appear- 
ing to arise from a perfect image of the source, in the 
three worst cases the apparent source is distorted and has 
moved toward the plane. In most practical situations, the 
agreement shown for the three best cases would be 


adequate, but it is significant that even at the furthest 
distance % = 160 A (7D 2 /A of the feed) there is noticeable 
error in the results. 

This error may be eliminated entirely by the use of a 
representation of the incident fields that is valid in the 
near field. Spherical waves are a well known set of 
solutions to Maxwell’s equations that satisfy this require- 
ment. If the incident fields satisfy Maxwell’s equations, 
then coefficients a n and b n may be found such that 

Hi (p, 0, <j>) = ~~r~ V a n n„ + b n m n (15) 

JOifi jLJ 
n - 0 

where n n is a transverse electric (TE) and m n is a trans- 
verse magnetic (TM) spherical wave, as given by Eq. 
(A-l) in Appendix A. 

The method for obtaining coefficients so that the wave 
expansion will match an arbitrary input pattern in the 
far field is also discussed in Appendix A; therefore, it will 
simply be assumed that coefficients have been obtained 



108 144 180 

POLAR ANGLE 0, d«g 

Fig. 6. Computed phase patterns for plane reflector 
cases; far-field assumption for incident fields 
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to describe the incident field. The form of Eq. (15) com- 
pletely eliminates the approximation of Eq. (14) by 
directly yielding H \ on the surface. However, the eval- 
uation of the spherical-wave expansion does involve 
computer time. Although this time is typically less 
than 10 percent of the total, the comment made earlier 
with regard to the integral-equation technique for ob- 
taining the currents also applies here: the expansion 
should be used only when the far-field form is inadequate. 

Clearly, the far-field form is inadequate in cases z — 5, 
10, and 20 A of the data presented in Figs. 4 and 6. These 
cases have been recomputed using a spherical-wave 
coefficient specification of the incident fields. The re- 
flected patterns are virtually indistinguishable from the 
incident pattern, as shown in Fig. 7. Phase data, trans- 
formed to the image point, are shown in Fig. 8; clearly, 
the error has been eliminated. 

Now that it has been established that the far-field 
assumption yields incorrect results that are corrected 
when a spherical-wave expansion is used, it is of interest 
to consider the physical mechanism of the discrepancy. 
Since the theory is based on current integration, the 
surface current distribution must reflect the essential 
difference between the two methods. The amplitude and 
phase of the y component (the principal polarization) of 
the surface currents at = 0 are compared for three 
cases in Figs. 9 and 10. The currents in Figs. 9 and 10 
are plotted versus the Cartesian coordinate x, rather than 
the polar angle 0 (see Figs. 1 and 3), to illustrate the 
true relative spatial distributions of the three cases. 

When the far-field approximation is assumed, the 
diameter of the region containing significant current mag- 
nitudes linearly approaches zero as the source approaches 
the reflector. Typically, a smaller source will result in a 
broader radiation pattern, but in this case the phase 
of the surface currents becomes more uniform as the 
source becomes smaller and this fact tends to make the 
pattern narrower. Over a considerable range of z values, 
these effects nearly cancel out, as demonstrated by the 
resulting radiation patterns shown in Fig. 4. However, 
in the limiting case of a point source, a broad dipole 
pattern would be expected with the phase center of the 
pattern on the reflecting surface. This is in fact exactly 
the trend exhibited by the data for the far-field case in 
Figs. 4 and 6. As mentioned previously, the surface cur- 
rents radiate symmetrically on both sides of the plane, 
and because a broader pattern cannot result in perfect 
destructive interference with the incident pattern, the 
result is the back lobe which appears in Fig. 4. 


When the spherical-wave expansion is used, the phase 
behavior is very similar to the far-field result, except that 
the 180-deg discontinuities are smoothed out. (Since the 
phase data are modulo 360 deg, these discontinuities 
could be in either the positive or negative direction. They 
were drawn to match the spherical-wave phase data as 
closely as possible.) The major difference between the 
spherical-wave and far-field cases is that the patterns of 
the current magnitudes are broadened, particularly in 
the z = 5 A case. Because the phase patterns are nearly 
the same, the larger source produces a narrower radiation 
pattern that brings the result into agreement with the 
incident pattern, as shown in Fig. 7. 

It is of interest to note that the radial field component 
of the incident magnetic field — which is completely ne- 
glected in the far-field approximation— makes a noticeable 
contribution to the total current, as shown in Fig. 11. 
When 2 = 5 A, this current contributes about 3 percent 
of the field strength of the scattered pattern on axis, so a 
small but not negligible part of the correction to the 
far-field approximation is the inclusion of the radial field 
components. However, these surface current data lead to 
the conclusion that the primary cause of the error result- 
ing from the use of the far-field approximation is that 
the current magnitude patterns are too narrow. 

The computer program developed in this report (see 
Appendix D) includes a subroutine FIELDS that pro- 
vides the main program with values of H* on the surface 
integration grid. Two such subroutines were written. 
In one subroutine the incident fields are specified by a 
set of spherical-wave coefficients a n and b n . In the other 
subroutine the far-field relations of Eq. (14) are assumed, 
and E 0 and E <p are represented by exactly the same type 
of series as the surface data (see Eq. 12). In the case of 
the field data, the m — 1 azimuthal variation term is 
of special significance and is the only term present in an 
important class of sources (Ref. 18). For this reason, the 
FIELD subroutines were written for any single value m 
for the order of azimuthal variation. If more than one 
value is required, the cases may be run individually and 
the results superimposed. The far-field subroutine is 
intended for cases in which the reflecting surface is 
sufficiently removed from the source, or in which the 
available number of spherical waves (limited by storage 
and other considerations) is inadequate to expand the 
incident-field pattern. The spherical- wave subroutine 
should be used whenever the reflector is close to the 
source, and in borderline cases (which include most 
practical Cassegrainian antenna systems) when assured 
accuracy is the overriding concern. 
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Fig. 7. Computed amplitude patterns for plane reflector cases; spherical-wave expansion of incident fields 


V. Evaluation of the Integral 

The numerical evaluation of the integral, Eq. (7), may 
justifiably be considered the most important part of the 
scattering problem. A primary engineering function of a 
scattering program is the evaluation of the effect of one 
or more parameters on the resulting pattern; this appli- 
cation is severely restricted if each case requires 6 or 8 h 
of computer time. 


in each of the variables of integration. In the case con- 
sidered here, the integral to be evaluated (see Eq. 7) may 
be written as follows : 


1(0,*) 


Mi 

2ir 





(16) 


where 


The fast Fourier transform technique has been applied 
to the evaluation of diffraction integrals by Pratt and 
Andrews (Ref. 19), but this formulation is applicable 
only when the argument of the exponential term is linear 


y = p * i r- p 

— p [sin 0 sin 0 (cos <j> cos ® + sin <j> sin <£) 

+ cos 0 cos 0 - 1] (17) 
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Fig. 8. Computed phase patterns for plane reflector 
cases; spherical-wave expansion of incident fields 


Since y is a nonlinear function of 6 and <j>, the fast Fourier 
transform technique is not applicable. Techniques that 
have been successfully applied to Eq. (16) are Gaussian 
quadrature, which was used by Allen (Ref. 20), and 
Romberg quadrature, which was recently used by Rusch 
and Strachman (Ref. 21). The Romberg method is par- 
ticularly useful in the case of computing the main-beam 
and near-sidelobe regions of a high-gain pattern, and 
could in fact be used with the method presented below 
to increase efficiency in this situation. 

Before describing the integration technique developed 
by the author, some general features involved in evaluat- 
ing this integral will be discussed. Suppose I(0,4>) is to be 
evaluated on a set of points (the output grid) given by 
(0 ; -, $*), 1 < / < JMAX, l<k< KMAX. Typical values 
are JMAX — 91, KMAX — 2, The integral will be eval- 
uated numerically when the integrand is specified at a 
set of points (the integration grid) given by (0 m> </> n ), 
1 < m < MMAX, 1 < n < NMAX. A modest grid size is 
MMAX = 50, NMAX = 181. (The effect of these grid 


parameters on computed pattern accuracy is discussed 
below and is considered in detail in Ref. 22.) 

As mentioned in Section II, with the far-field assump- 
tion, F does not depend on 0 and <5, so it need only be 
computed once on the integration grid — 9050 points in 
this example. This operation requires little computer time, 
but since F has three Cartesian components (each of which 
is complex), 6 X 9050 = 54,300 values must be stored. 
Many computers do not have this capacity. Therefore, the 
program was written so that several integration grids may 
be specified — i.e., the reflector is divided into segments — 
and the resulting scattered fields from each segment are 
superimposed. 


Unavoidably, the function must be computed at 
JMAX X KMAX X MMAX X NMAX points (1,647,100 in 
this example). This example demonstrates the tremendous 
advantage of storing F rather than recomputing it for 
every output point. Also, the number of points (1,647,100) 
represents the number of summations involved to evaluate 
the integral for all of the output points. Several steps were 
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Fig. 10. Surface currents ©n plane reflectors; 
phase behavior 



NORMALIZED COORDINATE, x/A 

Fig. 11. Surface currents on plane reflectors; contribution 
of the radial component of the incident fields 


taken to minimize the time used in this portion of the 
program: 

(1) The trigonometric functions sin 6 \ cos 6 , sin <£, cos </>, 
sin 0, cos 0, sin $, cos $ are precomputed at the 
required grid angles and stored so that the compu- 
tation of y will involve only multiplication and 
addition of stored quantities, rather than the time- 
consuming evaluation of trigonometric functions. 
For this reason, Eq. (17) is not written in the form 
cos </> cos 4- sin <£ sin <T> = cos (</> — $), which 
would have to be computed at every permutation 
of <f> and $ values. 

(2) A fast machine-language subroutine was written to 
compute cos ky + / sin ky. This subroutine is de- 
scribed in Appendix B. 

(3) A special numerical integration technique was 
developed that reduces the required values of 
MMAX and NMAX. 

The numerical integration technique is described in 
Ref. 5, but for completeness the description will be par- 
tially repeated here. For a fixed output point (® p , <f> q ) the 


integral is of the form 

f ®AI f 

/ / F(0 y <j>) exp jky PQ (#,</>) d0dcf> (18) 

J J fa 

where 

F = + psmeH^ji e 

+ p sin 6 He H P ) i $ 

+ (j^- sin ff„) i P (19a) 

H P i P + He i. + H^if = pe*P H» (19b) 

Since the unit vectors of Eq. (19) vary over the region of 
integration, Cartesian components are used in the actual 
evaluation of the integral. 

In the near field, the radial variation of H > will approxi- 
mate e~ ikp / i p behavior; therefore, if this dominant variation 
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is factored out, H p , H e , and H<p as defined in Eq. (19) will 
be nearly constant with respect to p, which makes it much 
easier to numerically approximate F. The e~ jkp factor has 
been included in the path length function y, and the 1/p 
factor is canceled by a p factor in the term for incremental 
area dS. Deviations from e~ jJip phase behavior are included 
by allowing H p> H 0? and H# to be complex valued. (Since 
a radial component of H* is included, this does not involve 
any assumptions or approximations). 

Consider the behavior of the integrand over an incre- 
mental area of S, 

AS mn — {($,</>) ' 0 m , (9 fC 0 m+ 1, <f> n <f) fC <pn+ 1 } (20) 

Suppose that the physical dimensions of AS m „ are on 
the order of a wavelength, A — 2 ir/k. Then the path length 
term jky cannot vary by more than 2tt, and since electro- 
magnetic fields cannot generally change abruptly over 


distances on the order of a wavelength, F will not vary 
much. Thus, in virtually any problem, F and y will be 
very well behaved and slowly varying over AS mn . How- 
ever, the possible 2tt variation in the exponential term 
could cause the real and imaginary parts of the integrand 
to behave like a full cycle of a sinusoid. To apply a tech- 
nique such as Simpson s rule to the entire integrand would 
require a further subdivision of AS. However, if the func- 
tions F and y are approximated individually a simple linear 
form will operate satisfactorily over AS m „. Explicitly, 
write 

^(^ 5 $) — &mn + bmn @m) ~"k (</> <f>n) (21a) 

— fynn “k fimn (0 $m) ~k £mn (<£ <£n) (21b) 

for (0 y (j))eAS mi . 


The method used for determining the coefficients is to best-fit a plane (least squares sense) to the values of 
the function at the corners of AS mn . For example, in the case of the function F this method yields the following: 


&mn ^ F (@ m +i,<f>n+i) ~k F(# m+1 ,<£ n ) “k F(0 m ,<£ n+ i)] 

-i 

[F(0m+i,<£») F($ w ,0 9 ) “I" F(0 H j + ij<£ n+ i) F(0„„*, 1+1 )] 

[F(0m,<t>n+i) -|- F(d,„ +1 ,<p n+1 ) F(<W»)1 


2A6>„, 

1_ 

2A fa 


(22a) 

(22b) 

(22c) 


where 


Aft* = ft* 


A <f>n — 




An identical relation holds between the coefficients a,,,,,, &»» and the function y. The integration over AS,, 

may then be performed analytically, yielding the contribution: 


aIjhh exp j l\CL mn 

A0, 


exp jkft mn A0 m 1 


+ b-m 


“k 


jkfimn 
exp jkfi nm A0 m 


. ikp 

exp jk/3 mn A0 m 1 
jkfimn 


- 1 " 

“ exp jkUA4>n ~ 

1 " 


- 

_ jk(z mn 

- 


/ exp jkf3 mn A0 m 1 \"1 

r exp jk£ m ,A<f>n — 11 

V 

(jkft mn ) 2 / J 

L 

jk£mn J 


A cf> n 
jk£mn 


exp jk£ m 


jk£ mn A(f>n — ^ 


exp jk£ mn A<f> n — 1 
(iUmn) 2 


(23) 


Since it is possible for /S mn or £ nm to be near or equal to zero, it is necessary to develop separate equations for 
this case (they are easily derived from the above equation) to avoid large numerical errors. 
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It should be noted that many of the computations 
involved in evaluating each single vector component 
do not have to be repeated for the remaining five com- 
ponents. For example, in Eq. (23) only the components 
of the ama, hmn and c mn coefficients will change. There- 
fore, although the six vector components present a severe 
storage problem, the effect on computer time is con- 
siderably less than a factor of six. 

The basic idea behind this method — to isolate the oscil- 
latory behavior of the integrand — is conceptually similar 
to the Eikonal technique used in deriving geometrical 
optics from Maxwell's equations (Ref. 23) and is also a 
feature of Filon's method for evaluating Fourier integrals 
(Ref. 24). In fact, Allen has also applied Filon's method 
to antenna problems and concluded that it is competitive 
with Gaussian quadrature (Ref. 20). The essential dif- 
ference of the method presented here is that the “fre- 
quency” terms J3 mn and are re-estimated for each 
incremental area of integration, rather than assumed as 
known constants. 

If radiation integrals are considered as analytic forms 
of a Huygens type of principle in which infinitesimal 
electric or magnetic dipoles radiate a simple pattern and 
are summed as in an array, then another way of inter- 
preting this technique is to describe it as replacing in- 
finitesimal elements with elements about a wavelength 
square that radiate a more complicated pattern as given 
by Eq. (23). 

Although this technique appears obvious and has an 
intuitive appeal from an engineering standpoint, it is 
somewhat unusual mathematically. The integral is a linear 
operator, and virtually all quadrature formulas are also 
linear, except this one. The easiest way to show this is to 
note that the functions 1 and e ike are both of the form for 
which the technique will be exact. However, their sum 

1 + gifce — 2 cos -g- e jke/2 (24) 

is not of the form (a + bO) e jk{0L+ & 6 \ and the sum of the 
(numerical) integrals is not equal to the (numerical) inte- 
gral of the sum. However, if F and y are well behaved, 
the technique does in fact converge to the integral oper- 
ator (and therefore becomes linear) in the limit of zero 
step size -(Ref. 5). 

Convergence tests for the case of scattering from a 
hyperboloid show that with this technique, incremental 
areas 2/3 of a square wavelength in size result in errors 


more than 40 dB below the pattern maxima (Ref. 22). 
Similar tests of a program that uses Simpsons rule to 
evaluate a one-dimensional integral show that the areas 
must be at most 0.04 square wavelengths for the two- 
dimensional case, and possibly smaller (Ref, 22). There- 
fore, the number of integration points may be reduced 
by at least a factor of 16 relative to Simpson's rule. The 
plane reflector data presented in a previous section were 
computed with the step sizes given in Ref. 22. These 
data are an absolute test of numerical accuracy and are 
further confirmation of the convergence tests. 

Although the complicated form of Eq. (23) seems to 
indicate that this technique would require an order of 
magnitude more time per data point than Simpson's rule, 
the fact is that in this type of problem the evaluation 
of e jx = cos x + j sin x tends to dominate machine time. 
For example, an IBM 7094 will multiply or divide two 
numbers in about 10 /as, or add two numbers in about 
15 fis (Ref. 25). However, computation of the sine and 
cosine to obtain the real and imaginary parts of e jx re- 
quires 571 /as with the library subroutine, or 281 /as with 
the fast subroutine described in Appendix B (see 
Table B-i). 

It is estimated that time per data point for this tech- 
nique is increased by a factor of 2 to 4, relative to 
Simpson's rule, which means that the net reduction in 
total running time is a factor of 4 to 8. In practice, this 
technique can result in a reduction from 8 h to 1 or 2 h 
of computer time, so it is of substantial practical impor- 
tance. For example, Slobin (Ref. 16) has reported that a 
program using Simpson's rule integration requires 1.3 min 
of IBM 360/75 computer time per output point to com- 
pute the scattered pattern of a 56-wavelength-diam tilted 
subreflector. Thus, 3.90 h would be required for 180 
output points. These data have been obtained in 1.01 h 
of IBM 7094 Mod I computer time with the program 
described here. (The 1.01 h includes 9.8 min to evaluate 
the spherical-wave expansion of the incident fields at 
each of the integration grid points.) The 360/75 is faster 
than the 7094 Mod I by a factor of between 3 and 5; 
therefore the program described here is 11 to 19 times 
faster than this particular program that uses Simpson's 
rule. (It is interesting to note that the new technique 
works better on an old computer than the old technique 
on a new computer.) 

It should be noted that in the case of rotational sym- 
metry, when one integration may be performed analyti- 
cally, these same data can be obtained by evaluating a 
one-dimensional integral (with a Bessel function in the 
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integrand) in about 6 min (Ref. 22). Therefore, it is 
grossly inefficient to evaluate a double integral in sym- 
metrical cases, and the program described here should 
only be applied to problems involving asymmetrical con- 
figurations. 

These convergence data are representative of the “low- 
gain” case in which the contributions from currents on 
the surface do not add in phase at any output point. A 
special case of considerable practical interest is the com- 
putation of the main lobe and first few sidelobes of the 
scattered pattern of a high-gain antenna such as a para- 
boloidal reflector. In this case, contributions are adding 
nearly in phase, i.e., the phase of the integrand is venj 
slowly varying. Tests recently completed indicate that, in 
this case, integration points may be spaced at intervals 
of 20 wavelengths or more without seriously affecting the 
accuracy of the results. (Actually, in many practical cases 
closer spacing will be required to represent sufficient 
detail of the surface and incident-field data.) The results 
of these recent tests are similar to results reported by 
Allen for Filons method and Gaussian quadrature for 
this case (Ref. 20), so the integration technique used here 
is nearly equivalent to the techniques in this situation. As 
mentioned previously, any of these methods may be im- 
proved by a Romberg type of modification in this case. 
The large spacings allowable between grid points means 
that the patterns of reflectors over 1000 wavelengths in 
diameter may be computed. However, beyond the first 
few sidelobes the situation begins to approach the “low- 
gain” case described earlier. 

VI. Comparison With Experimental Data 
and Other Results 

An improved theory is important if it confirms the 
results of an approximate theory, since the approximation 
may then be used with confidence; it is also important 
if an improved theory corrects the results of an approxi- 
mate theory. The application of spherical-wave theory 
developed in this report is considered valuable on both 
counts. A case in which the far-field approximation yields 
an incorrect result has already been presented; a case in 
which this approximation is adequate, except for the 
most precise calculations, is discussed in this section of 
the report. 

A basic criterion in selecting this sample case was that 
it represent a real problem. The standard NASA/JPL 
Deep Space Network antenna is an 85-ft-diam parabo- 
loid, with a Cassegrainian feed system, operated at 


S-band. This configuration is also common for communi- 
cations satellite terminals, so this sample case is clearly 
one that is of practical interest. 

The object of the analysis was to compute the radiation 
pattern of the overall antenna, with particular emphasis 
on the pattern maximum that determines the antenna 
gain. The antenna is illustrated in Fig. 12. The design of 
this particular antenna is well documented (Ref. 26), and 
JPL provided experimental patterns of the fields scat- 
tered by the subreflector used with this configuration. 

There is reason to question the application of the far- 
field approximation to the analysis of this antenna, because 
the subreflector is at 1.65 D 2 / A of the primary feed, and 
the main reflector is at 0.148 D 2 /A of the subreflector. 
Because the results of a far-field analysis generally agreed 
with experimental data, it was certain that the approxi- 
mation was reasonably good, but the amount of error 
introduced was uncertain until the completion of the 
analysis presented in this report. 
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The primary feed for this antenna was the same dual 
mode conical horn (Ref. 27) used for the plane reflector 
cases. The spherical- wave coefficients for this feed were 
obtained from an experimental pattern using the tech- 
nique discussed in Appendix A; the incident-field pat- 
terns shown in Figs. 7 and 8 represent an evaluation of 
the spherical-wave expansion at p = oo . Two parallel 
analyses were performed on the antenna. One analysis 
was made with the spherical-wave expansion FIELDS 
subroutine and the other with the far-field approximation 
FIELDS subroutine. The analyses were identical in all 
other respects. 

The two computed scattered patterns of the subre- 
flector (which consists of a vertex plate^ a hyperboloidal 
section, and a conical flange) were virtually identical. 
This is not surprising, because the subreflector is near 
the 2D 2 /A distance from the primary feed. 

The experimental subreflector pattern was measured 
at a range that was nearly equal to the distance to the 
main reflector. In Fig. 13, these experimental data show 
good agreement with the pattern computed with the far- 
field approximation. The pattern computed using a 
spherical-wave expansion was in turn expanded in a 
series of spherical waves. In Fig. 14, the spherical- wave 
expansion, evaluated at the same radius at which the 
experimental data were measured, is compared to the ex- 
perimental pattern. The agreement is even better than 
that in Fig. 13. The computed data closely follows the 
experimental pattern, demonstrating that the agreement 
between theory and experiment is outstanding. 

The next step in the analysis was to use the subreflector 
pattern to illuminate the main paraboloidal reflector, 
again by the use of a spherical-wave expansion and the 
far-field approximation in the respective cases. The two 
computed main reflector patterns are compared in Fig. 15. 
Although the most noticeable difference is in the side- 
lobes, this difference will generally be of little impor- 
tance. Of greatest significance is that the peak of the 
main beam (i.e., the computed gain) differs by 0.127 dB 
in the two cases. 

With the scale used in Fig. 15, a 0.127-dB difference 
is difficult to detect, and in many applications of com- 
puted antenna data it will be of no importance whatso- 
ever. Therefore, for these applications (and for antennas 
of this approximate design) it has been shown that the 
faster and simpler far-field analysis is inadequate. For 


example, the far-field analysis will be used for the feed 
displacement calculations to be presented shortly. Use 
of the far-field analysis is of particular importance in the 
case of development work involving a large amount of 
computer time. 

However, there are cases in which a 0.127-dB differ- 
ence is significant. Computed data have reached a level 
of sophistication where they are frequently weighted 
about equally with experimental data in calibrating large 
ground antennas (Ref. 28). An 85-ft antenna costs about 
one million dollars (Ref. 29), and it can be argued 
(Ref. 30) that an uncertainty of 0.127 dB in the gain of 
such an antenna is worth $40,000. In fact, the costs of 
the time and effort spent to calibrate an antenna with 
this precision are even greater than $40,000. Since the 
IBM 7094 computer time required for the two cases 
shown in Fig. 15 (including the spherical- wave expansion 
of the subreflector pattern) was 23.24 and 8.60 min for 
the spherical-wave and far-field analysis, respectively (a 
difference of about $61.00 at the current price of com- 
puter time) it is probably safe to say that most engineers 
would choose the more accurate analysis. (Since the eval- 
uation of the spherical- wave expansion on the integration 
grid is independent of the number of output points, the 
relative times would be more nearly equal for a larger 
number of output points.) 

In all of the cases evaluated so far, rotational sym- 
metry has been assumed. As a last example, the problem 
of an offset feed system will be considered to illustrate 
some of the capabilities of the program for asymmetric 
geometries. 

The case considered is the offset feed geometry shown 
in Fig. 16. To provide a basis of comparison, the dis- 
placements along the dashed curve shown in the figure 
are the same as a case previously analyzed by Ruze 
(Ref. 31). However, the true subreflector pattern has 
been retained, thereby eliminating the need to assume 
a simple form for the illumination of the main reflector. 
In addition to the displacement translations, the feed 
system was rotated so that the edge angle of the main 
reflector remained nearly constant, to minimize spillover. 
The offset angle is measured in the conventional 
manner. The computed patterns for a zero offset and 
two nonzero offsets are shown in Fig. 17. The analysis 
given by Ruze, which is based on scalar theory, includes 
the case of an illumination function of the form 
f(x) = 0.3 + 0.7(1 — x 2 ), that approximates the reflector 
illumination in the case considered here. 
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POLAR ANGLE ©, deg 

Fig. 13. Comparison of experimental and computed subreflector patterns; far-field approximation case 


For this geometry and illumination, Ruze predicted However, in addition to the confidence factor mentioned 

a beam deviation factor of 0.84. The beam deviation earlier, the approximate analysis has several restrictions: 

factors of the results given in Fig. 17 are 0.835 and 0.840, (1) it assumes a simple analytic illumination function, 

for the 1- and 2-deg cases, respectively. The coma lobes (2) the displacements must be along the specific curve 

agree within 1 dB with the values predicted by Ruze and shown in Fig. 16, and (3) the displacements are assumed 

the half and tenth-power beamwidths agree within to be small compared to the focal length. In the case 

0.04 deg. The gain loss agrees (as well as it is possible to considered in this section, these restrictions were met by 

read the graphs given by Ruze) within approximately design to obtain a comparison with independent data, 

0.1 dB. but the program is capable of a case with (1) a very 

complicated illumination pattern (for example, a pattern 
with a central null produced by a vertex plate, including 
This agreement is so complete that it might seem the phase pattern perturbations typical in this situation), 

irrelevant to bother with the more rigorous analysis. (2) the displacements in any direction and, (3) a tilt 
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POLAR ANGLE deg 


Fig. 14. Comparison of experimental and computed subreflector patterns; spherical-wave expansion case 


angle that is restricted to approximately 5 deg in the 
present program (the wavelength determines the allow- 
able peak error that in turn determines the tilt angle as 
shown in Fig. 2), but that could be made arbitrarily large 
if more Fourier components were included, or if a spe- 
cial surface subroutine were written for the particular 
case of a tilted paraboloid. 

VII. Summary 

The results shown in the previous section clearly 
demonstrate the usefulness of the physical-optics tech- 
nique. The agreement between theoretical and experi- 


mental data is excellent, and in the case of the 85-ft 
antenna, the application of a technique such as the 
integral-equation method to obtain the surface currents 
is virtually hopeless because of the large size of the 
reflector. 

It was demonstrated in Section IV that the far-field 
approximation for the incident fields can lead to poor 
results, and that the spherical-wave representation cor- 
rects the results. The far-field expansion is still useful in 
many practical cases, but if precise results are required, 
and if the reflectors are in or close to the near-field 
region, the spherical- wave representation should be used. 
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coincident with a coordinate boundary, an approximate 
boundary must be constructed out of trapezoidal-like 
segments. An interesting extension of the integration 
method developed here would be to compute the scat- 
tered pattern of a triangular-shaped segment to be used 
in constructing edges that are at least continuous. An- 
other possible improvement is to modify the program to 
accept an arbitrary number of Fourier components in the 
surface representation that would eliminate the small 
restriction on the reflector tilt angles. 


The main area of future work lies in applying the 
scattering program. Parameter studies of near-field 
Cassegrainian antennas using actual experimental feed 
patterns is an obvious example. In all of the cases con- 
sidered above, the spherical-wave coefficients of an 
incident-field pattern were obtained with the program 
described in Appendix D, but the scattering program will 
accept coefficients regardless of how they are obtained. 
For example, if the currents excited on a feed device, 
such as log-periodic wire antenna, are known (perhaps 
from the integral equation method), the spherical-wave 
coefficients may be obtained directly from the currents. 


Fig. 15. Computed patterns of 85-ft-diam main reflector; 
comparison of far-field and spherical-wave expansion 
cases 


A Fourier type of representation of the scattering sur- 
face is efficient in the case of a slightly tilted figure of 
revolution, and in many cases three components are 
sufficient to accurately specify the surface. 

In the case of a low-gain asymmetrical reflector the 
program developed here is an order of magnitude faster 
than an earlier program; in the high-gain case the inte- 
gration technique loses much of its advantage over prior 
methods, but the other timesaving techniques used still 
result in efficient operation. For symmetrical geometries, 
one integration should be performed analytically. 

VIII. Possible Improvements and Future 
Applications of Computer Program 

The computer program could be improved in several 
ways. For example, a Romberg type of procedure could 
be incorporated into the integration method. Also, in the 
existing program, when the edge of the reflector is not 
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Fig. 17. Computed patterns of 85-ft-diam main 
reflector with offset feed 


Also spherical-wave coefficients have been obtained 
analytically for cases such as a plane wave and an arbi- 
trarily located and oriented dipole (Ref. 32). This dipole 
case is precisely the expansion required to compute the 
fields in the focal region of a reflector by the use of 
the same reciprocity technique recently used by Rusch 
(Ref. 21). This is a very interesting area for future work, 
since focal-region fields studied have involved assump- 
tions and restrictions, whereas this method is, in principle, 
exact. Also, this method is applicable to nonparaboloidal 
antennas, such as shaped Cassegrainian antennas 
(Ref. 33), that have never been analyzed to obtain the 
focal-region fields. The reciprocal part of this problem is 
the computation of the pattern of reflectors with de- 
focus ed feeds, and this is also a very interesting area for 
application of the scattering program, as discussed in the 
previous section. 

The spherical-wave expansion program also has future 
applications independent of the scattering program. For 
example, a study is currently being pursued to apply 
spherical-wave expansions to the determination of cor- 
rection factors for near-field gain measurements (Ref. 34). 
The usefulness of mode expansions in the solution of 
boundary-value problems has not been exhausted, and 
cases may be found in which the method is superior to 
the currently popular integral-equation method. 
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Appendix A 

Spherical-Wave Expansions of Electromagnetic Fields 


In a source-free, homogeneous, isotropic medium, 
spherical-wave solutions of Maxwell’s equations are given 
by Stratton (Ref. 35): 


q= z n (kp) 


mP n (cos 6) sin 


sin# 


cos 


m<t> i e 


- *» ( k p) p n (cos 6) m<f> (A-la) 


n* mn — 


n(n + 1) P m ( cos d ) sm m<f> ^ 


+ T^l^ [ P z ^ k P )] ie P " icos6) cos m4>i ° 


-LJ-[p Zn (kp)]^^ 


mP »‘ (cos e ) cos ,/n 


Sill 


771(f) l<p 

(A-lb) 


where 

e jtat = time dependence (implicitly 
assumed) 

(p, #,</>) — spherical coordinates (see Fig. 1) 

k = o)(e/x) 1/2 = 2 tt/A = propagation constant 

% n (kp) = any solution of the spherical Bessel 
equation 

P" 1 (cos 6) ~ associated Legendre function 


In the following derivation, an electromagnetic field is 
understood to be a pair of vector valued functions E 
and H (defined everywhere in a source-free, linear, ho- 
mogeneous, isotropic region V) that satisfies 



(A-2b) 


are defined and finite. 

Since V • V X A = 0 for any vector A, Eq. (A-2a) is 
equivalent to the assumption that Maxwell’s equations 
are satisfied; by Poynting’s theorem, Eq. (A-2b) is equiva- 
lent to the requirement that sources radiate finite power 
(assuming finite energy in any resonant fields). 


The parallels between electromagnetic theory and the 
theory of complex variables are intriguing, and it is 
interesting to note that by Eq. (A-2a), the existence of the 
first derivative ( V X E and V X H) assures the existence 
of derivatives of all orders (VXVX---XVXE, etc.). 
Other striking analogies are between the field integrals 
(which express the fields everywhere in a volume in terms 
of the values on the surface enclosing the volume) and 
Cauchy’s theorem, and between spherical-wave expan- 
sions and Laurent series. 


The objective of this section is as follows: given an 
electromagnetic field in a region V that consists of all 
space outside of a sphere of radius p 0 (i.e., all sources are 
enclosed in this sphere), determine coefficients a e omn and 
be mn so that everywhere in V 

'y y y j &Q?nn 7K\ e 0 mn "L b^ mn tl^ mn (A-3a) 


V X E and V X H exist everywhere in V 


H(pM) = 7-7Z Z a lmn n Inn + 
f 0 i r J m n 




(A-3b) 


and 


V X E — — 


V X H = /WE 


(A-2a) 


where e and /a are physical constants, and to is a nonzero 
positive constant, and where 


Jones (Ref. 32) has shown that any electromagnetic 
field can be written in this form. Furthermore, since 
sources were assumed to be of a finite extent, only the 
solutions involving h™ (kp) (the spherical Hankel func- 
tion of the second kind that satisfies the radiation condi- 
tion and corresponds to an outward traveling wave, as 
discussed in Refs. 36 and 37) need be included. 
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Bouwkamp and Casimir (Ref. 38) obtained expressions 
for the coefficients in terms of the source currents; 
Kennaugh (Ref. 36) gives expressions in terms of the 
tangential E and H fields on a surface enclosing the 
sources; Jones (Ref. 32) includes expressions in terms of 
all components of either E or H. Since the data involved 
here are the tangential components of E on a sphere of 
radius p t > p 0 , a derivation for this case will be given. 
(Note that the case Jones considers is equivalent for the 
determination of TE wave coefficients but not for the de- 
termination of TM wave coefficients.) 

As a first step, let a^ mn and b e Qinn be any set of coeffi- 
cients so that Eq. (A-3) is satisfied. (As mentioned above, 
Jones has shown at least one such set must exist.) 

Define 


0*m» = %mn K* ( kp 1 ) (A-4a) 

b\ mn = (kp)] P=Pl (A-4b) 


Then, since Eq. (A-3) is true everywhere in V, in par- 
ticular, there must be equality of the tangential com- 
ponents at p = p t : 


Eq ( 0,<f > ) ~+~ 


mP™ (cos 0 ) S in 
cos m <t> 


sin0 


+ b tmn -^P>:(cose) 


(A-5a) 


E*(e,<t>) =££ 


a 'omnjj P: (COS 0)';r:m<l> 


cos , 
sin 


m ?;; 1 (cos 0 ) cos 


sin 0 


sm 


m<t> 


(A-5b) 


where 


E e = E (p u 0,<f>) • i q 
E,p(e, 4 >) = E( P M ' U 


The azimuthal components are separated out by the 
use of an ordinary Fourier expansion. (These equations 
are true for m > 0. For m — 0, an additional factor of 2 
should be included but has been omitted for clarity.) 


r 27T 

A< ^ “ Jo Ee ^ cos m4> d * 

„ ' mP“ (cos 6) 

= ~ * E =p a i ™ 


sin# 


+ b‘ 


9# 


(cos 9) 


(A- 6 a) 


T 27T 

(0) = / Ef (#,<*>) ^ m<f> dj> 


~ ae o™ 0? p » ( cos 9 ) 

+ c ™ p :: (cos e) 

~ ° mn sin 9 " 


(A- 6 b) 


For the next step, the following relationship is needed: 
mP!" (cos 9) 0 If mP"' (cos 9) 

— V, ±4r p, ;‘(cos 9) — v . — 

sm# 9# lV 'J L sin# 

1 . PT 9 mV 


i:i 

i'R 


m 2 P”'P” 


1 ~ » _9_ p I»_9_ p 

sin 2 # T 9# 1 9# ‘ 


-^•P“(cos #) si 


sin 6 dS ~ 


sin (9 d0 : 


— P w 
90 i sin 0 


+ 


mP? 


sin 0 00 


si 


P™ sin 6 dd — 


0 for l^n 

2 (n + m)! , 1X 

2 n + 1 (» - m) ! n ^ n + ^ 

The equation above may be obtained by the performance of one integration by parts 


n 


9 mP„ 

0 pm n 

00 7 sin 0 


+ 


mV" 


sin# 9# 


I 

n J 


sin # d# 


for i 


n 


(A-7) 


n o mP”‘l fT a mP™ 1 

sin * * + [mP T p “ ] o - 1 Lw p ”‘ inryj sin = 0 

and application of a result given by Stratton (Ref. 35) (Section 7.13, Chapter VII). 


(A- 8 ) 
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Then it is easily verified that 


jT[ ± m? ^ C ° S - ■ A m e o W + -^r P ? ( eos e ) B < (*)] sin ede = ” a W 21 + l (t-mji W + (A ' 9a) 


J o --£0 p T (cos e ) A m e (e) q= mPf ^ cos ^ 


sin 0 


. B m e (0)] sin edO = Tr U nl - jftgjj- i(i + 1) (A-9b) 


Equations (A-4) through (A-9) may be rewritten in the expected form: 


a *mt — 


1 2n + 1 (n — m) ! 

[Zn (fcpi)] 2 ir2n (n + 1) (H + m ) ! 


X 


n Tt 
» 


~*R e 0 mn * Ei(pi,0,<£) tangent iai sin 0 d0 d<f> (A- 10a) 


= 


2n + 1 (n — m ) ! 




?r2n (n + 1) (n + m) 


| /'21T fTT 

r x I /. 


“ n ®wm * Ei (/M, ^tangential sin (9 d0 dcf> (A-10b) 


Equation (A-10) not only gives an explicit form for the 
desired coefficients, but, by the well known properties of 
trigonometric and Legendre function expansions, the co- 
efficients are uniquely determined. Therefore, this illus- 
trates the fact that an electromagnetic field as defined 
earlier is uniquely determined by the value of the tan- 
gential E field on any sphere of radius pi > p 0 . 

The spherical-wave expansion program (orthogonality 
version) given in Appendix D evaluates Eq. (A-10) to 
obtain a^ mn and Z?* mn . That is, it is assumed that the 
azimuthal Fourier expansion of Eq. (A-6) has already 
been performed, and the input data to the spherical-wave 
expansion program are the (tabular) functions A m e (0) and 
B m e ( 0 ). Because successive cases with different values of 
m may be evaluated, in principle an arbitrary pattern 
can be expanded. However, the m ~ 1 case is of particular 
importance, as mentioned previously, which is why the 
program was written to handle only a single azimuthal 
component at one time. The odd components have also 
been neglected to avoid unnecessary complication, since 
the vast majority of problems may be analyzed using 
linear polarization. In the case of far-field input data, 


Eq. (A-4) may be rewritten with the asymptotic form of 
the Hankel functions, thereby yielding for p t oo 

= (A-lla) 

= (A-Hb) 

The scattering program implicitly assumes Eq. (A-ll) 
so it accepts the wave coefficients as output from the 
spherical-wave expansion program (it is assumed a 
far-field pattern is input to the spherical-wave expansion 
program). 

Since the solution is known to exist and because it is 
unique, it is also possible to solve Eq. (A-6) directly. It is 
well known that if all sources are enclosed in the sphere 
of radius p 0 , only modes for which n < k 0 will make a sig- 
nificant contribution to the field (Refs. 39-41). Therefore, 
Eq, (A-6) can be solved as a system of linear equations 
to obtain the coefficients. Again, with only even terms 
and a single m component, the m and ® subscripts can be 
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dropped and the matrix equation can be written: 


-A&n 


~ F“ (8 1 ) — G“ (8,) F»‘ (8 1 ) ••• -G» (8 1 ) “ 
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Qm fa) Qm (^) . . . ( 8 ,) 
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(A- 12) 


where 


F": (eh 

G?(e)i 

N t 


mP™ (cos 0) 
sin 6 


dO 

kp 0 


P™ (cos 0) 


Since (0 deg) = G” l (0 deg) and F™( 180 deg) — — G™ 
(180 deg), the values 0 — 0 and 180 deg would make the 
above matrix singular; therefore, they must not be used 
as data points. Even excluding these values, the first 
attempts to invert this system of equations indicated that 
the matrix was extremely ill-conditioned (on physical 
grounds it can not be singular if the 0 values are distinct 
and if 0 = 0 and 180 deg are excluded). It was sub- 
sequently determined that the matrix was ill-conditioned 
because the data points were concentrated in the region 
between 0 and 40 deg (where A(0) and B(0) had sig- 
nificant amplitudes). When the data points were equally 
spaced between 0 and 180 deg, the matrix was easily 
inverted. The program that was used to accomplish 
this, the spherical-wave expansion program (linear equa- 
tion version), is also presented in Appendix D. When 
properly used, both programs yielded the same results 
except for modes containing a negligible fraction of the 
total power. 


The linear equation method is included here primarily 
because it may be easily generalized to be valid when 
the fields are known on any surface, whereas the formu- 
lation of the orthogonality technique is valid only over 
the surface of a sphere. Furthermore, the linear equation 
method provides an interesting numerical check on the 
other program. However, the orthogonality method is 
more convenient for a number of reasons: (1) it is faster 
to perform the numerical integration than to invert the 
matrix, (2) the number of mode coefficients is not re- 
stricted to the number of input 0 values, and (3) Parseval’s 
formula can be used to show that a sufficient number of 
modes have been considered, rather than the well docu- 
mented but essentially qualitative arguments leading to 
the N ~kp relationship. 

The equations for obtaining wave coefficients derived 
here differ slightly from those obtained previously. The 
main difference is that these equations have been pro- 
grammed and used to obtain wave expansions of three 
kinds of patterns — analytical, numerically computed, and 
experimental. Previous work has been restricted to inter- 
esting but idealized cases in which the coefficients were 
obtained analytically in closed form; e.g., Potter (Ref. 41) 
found the coefficients for a circularly symmetric optimum 
illumination pattern; Kennaugh and Ott (Ref. 42) found 
the wave excited by a plane wave normally incident on 
a paraboloid; Jones (Ref. 32) gives the coefficients for a 
plane wave and for an arbitrarily located electric dipole. 
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Appendix B 

Fast Trigonometric Subroutine 


The machine-language function subroutine* EXPJX 
was written for the IBM 7094 computer and returns the 
complex-valued result e ix = cos x + / sin x , for real x , in 
one-half the time required by the IBM library subroutine, 
and with only slight reduction in accuracy. 

A table lookup technique is employed where 64 values 
of sin x m at x m — m (tt/128) + (tt/256); m — 0, 1, • • 63 are 
built into the subroutine. 

The basic algorithm used is as follows: 

(1) sin ( — x) = —sin ( x ), cos (— x) — cos (x), assume 
x > 0. 

(2) Write x = n (tt/128) + r (tt/128) where n is an 
integer, and 1 > r > 0. 

(3) Reduce n mod 256 to discard multiples of 2 ?r. 

(4) Write n — 64 <7 +* m, where 63 > m > 0. Then q + 1 
is the quadrant of the argument, and m is the in- 
dex of the tabular value x m closest to x (reduced 
to the 1st quadrant). 

(5) Define Ax = x — x m , = (r — 1/2) tt/128; and | Ax | 
< tt/265. 


(7) Use the quadrant index q and the sign of the orig- 
inal argument to obtain sin x and cos x for the 
unreduced (true) argument, and return the com- 
plex number (cos x, sin x). 

The reason for choosing 64 (or in general 2 n ) tabular 
values is that a single multiplication of the input argu- 
ment x( 128 /tt) results in a binary number in which the 
fractional part is r, the next six higher order bits are m, 
the next two bits are q , and higher order bits represent 
multiples of 2rr. Therefore, steps (2), (3), and (4) may be 
accomplished very quickly by one multiplication followed 
by shifting and masking operations. 

The MAP subroutine, which was named EXPJX, re- 
quires a total of 167 (decimal) storage locations, com- 
pared to 123 locations for the IBM Fortran IV library 
routine FCSN. 

An accuracy and timing test was made in which re- 
sults from EXPJX and FCSN were compared to the IBM 
double precision Fortran IV library routine FDSC. The 
library routine was called twice to obtain values for both 
sin x and cos x. Values were computed for 50 X 10 3 
random arguments uniformly distributed over the range 
(—5 to 5 tt), with results as shown in Table B-l. 


(6) Use cos x m — sin x 63 _ m , the expansions 

sin (x m 4* Ax) — sin x m cos Ax + cos x m sin Ax 

/ , Ax 2 \ / Ax 3 

^ smXtt ( 1 — -gf 1 + COS Xm ( Ax — -gj- 

cos (x m + Ax) — cos x m cos Ax — sin x m sin Ax 

/ 1 Ax 2 \ . / Ax 3 

— COS X m f 1 - -g|-J “ sm Xm f Ax gj" 


The relative accuracy of EXPJX becomes poor for 
small values of sin x or cos x, and for arguments | x | > Sir 
the loss of significance in the reduced argument will in- 
crease the errors, 

The subroutine of Printout B-l has been submitted for 
distribution through SHARE (SHARE library distribu- 
tion number SDA 3534), and COSMIC (Program NPO 
10439). 


and the approximation 

* “ -fn — ^ ~ 0,16 x 10 " 4 ) 

for | Ax | < (tt/256) to compute sin x and cos x in 
the 1st quadrant. 


*J. Hatfield of JPL Section 314, and H. Thacher of Notre Dame 
University were associated with the author in the development 
of this subroutine. 


Table B-l . Performance statistics for fast 
trigonometric subroutine 


Program 

Absolute 
error (rms) 

Maximum 
absolute error 

Average 

speed, 

MS 

EXPJX 

8.2 X 10“ 3 

3.7 X 10" 7 * * 

281 

FCSN 

0.28 X 10~ s 

0.13 X 10' 7 

571 a 

a For both sine and cosine. 
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Printout B-l. Computer printout of fast trigonometric subroutine 


$IBMaP expjx 
LBL 
TTL 
REM 
REM 
REM 
REM 
REM 
SPACE 
ENTRY 

N SET 

SPACE 

REM 

SPACE 

REM 

SPACE 

REM 

REM 

REM 

SPACE 

REM 

SPACE 

REM 

REM 

REM 

REM 

REM 

REM 

REM 

REM 

REM 

REM 

REM 

REM 

REM 

REM 


LIST, REF, M94, ( ) OK 

expjx, begin, 

COMPLEX EXPONENTIAL SUBROUTINE,. .3EXPJXS 

# i'f. % £ £ sjc # >Jt # # sjc jjc * # >|c # jjc # # ^ jjc # £ # # sji # # ajc s{c sje # $ # # # jJs# # # sJ»J< # * £ 

******* COMPLEX EXPONENTIAL SUBROUT I NE ... 3bXPJXa******* 
** * * * * * * * * ** * * 

* * * * * **** * * * * * * * * * * * * * * * * * * * * * * * * £ £ ^ ^ * * * * * * * * * * * sjc * * * * * * 

1 

EXPJX 

64 SIZE OF SINE-COSINE TABLE 

2 

THIS SUBROUTINE COMPUTES 
1 

EXP{ J*X)=COS(X)+J*SIN(X) 

1 

GIVEN THE FLOATING PT. VALUE, X. 

THE COS(X) AND SIN(X) ARE TABULAR EVALUATED AND SEXPJX3 
RETURNS WITH THEIR RESPECTIVE RESULTS IN THE AC. AND MO. 
3 

PROCEDURE... 

1. COMPUTE Y = ( 1/DTHETA ) #X AND SAVE SIGN UF X 

2. UNFLOAT Y WITH INTEGER IN AC. AND FRACTION IN MQ 

3. UNPACK INTEGER SO THAT QUAD. NO, AND TABLE INDEX 

ARE BITS 28-29 AND 30-36 RESPECTIVELY 

4. COMPUTE DX=ABS( FRACTION OF Y)*DTHETA-UTHETA/2 

5. COMPUTE T1=SIN(DX) AND T2=C0S(DX) USING SERIES 

EVALUATION (AT MOST TWO TERMS) 

6. COMPUTE Fl=SIN(XS)*Tl+COS(Xa )*T2 

F2=C0S( Xa))*Tl-SIN(Xb))*T2 
WHERE SIN(XS). AND C0S(X3) ARE OBTAINED FROM 
IT3S CORRESPONDING TABLE X3=X-DX 

7. FOR THE APPROPRIATE QUADRANT COMPUTE.. 

1ST QUAD. SIN<X)=(SIGN) FI 

COS ( X ) =F2 


REM 

REM 

REM 

REM 

REM 

REM 

REM 

SPACE 

EXPJX SXA 
SXA 
SXA 
SPACE 
LOG* 
STQ 
FMP 
SSP 
UFA 
RQL 
LRS 
ANA 
LGR 
PAC 
PXA 
LGL 
PAC 
PAX 
MPY 
SUB 
STO 
ARS 
CHS 
XC A 
MPY 
ALS 


2Nq 

quad* 

SIN(X) 

= C S IGN ) 

F2 



COS(X) 

= -F 1 


3RD 

QUAD. 

SIN(X) 

=- ( SIGN ) 

FI 



COS(X) 

=-F2 


4TH 

QUAD. 

SIN(X) 

=-( SIGN) 

F2 



COStX) 

= F1 



8. FLOAT SIN(X) AND COStX) 

3 

E X I T 1 , 4 
EXIT1 + 1, 1 
EXITl+2,2 
1 

3,4 FETCH X (IN RADIANS ) 

SIGN 

10VDT COMPUTE X/DTHETA 

=0233000000000 UNFLOAT IT (I IN ACC., F IN MQ.) 

8 SCALE F AT BIT 0 IN MQ. 

0 SET SIGN OF F POSITIVE 

=0377 MASK I FOR QUADRANT AND INDEX NUMBERS 

6 

0,1 XRl=QUADRANT number 

0,0 

6 

0,2 XR2= I NDEX NUMBER FOR SINE 

0,4 XR4= I NDEX NUMBER FOR COSINE 

DT COMPUTE FS) = F*(PI/2)/64 

DT0V2 

DX SAVE DELTAX=Fh)-(PI/2)/l28 

1 


DX 

1 CHANGE SCALE FROM BIT 2 TO BIT 1 
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Printout B-1 (contd) 



ADD 

= 1 B1 



STO 

T 1 

Tl=l-DELTAX**2/2 (COS( DE(_TAX ) ) 


LDQ 

$ I NX » 2 



MPY 

T 1 



STO 

FI 

F 1 = S I N { Xa ) *T1 


LDQ 

C0SX,4 



MPY 

T2 



ADD 

FI 



STO 

FI 

Fl=SIN(Xa)*Tl+C 0 S(Xa)*T 2 


LDQ 

S I NX * 2 



MPY 

T2 



STO 

F2 

F2 = S I N ( XS ) *T2 


LDQ 

CQSX ♦ 4 



MPY 

T 1 



SUB 

F2 



STO 

F2 

F2=COS(Xa)*Tl-SIN(Xa)*T2 


TRA 

*+1,1 



tra 

QUADl 



TRA 

QUAD2 



TRA 

QUAD3 



TRA 

QUAD4 



SPACE 

1 


QUADl 

CLA 

FI 



LDQ 

SIGN 



LLS 

0 



LDQ 

F2 



TRA 

EXIT 

(SIN(X)sSIGN.Fl, COS ( X ) =F2 ) 


SPACE 

1 


QUAD2 

LDQ 

F2 



CLA 

SIGN 



LRS 

0 



CLS 

FI 



XC A 




TRA 

EXIT 

( S I N ( X ) =S IGN* F2 , COS ( X ) =-Fl ) 


SPACE 

1 


QUAD3 

LDQ 

FI 



cl a 

SIGN 



CHS 




LRS 

0 



CLS 

F2 



XCA 




TRA 

EXIT 

{SiN(X)*-SIGN.Fl, COS ( X ) =-F2 ) 


SPACE 

1 


QUAD4 

CLA 

F2 



LDQ 

SIGN 



LLS 

0 



CHS 

LDQ 

SPACE 

FI 

(SIN(X)=-SIGN.F2, COS ( X ) =F 1 ) 

EXIT 

STO 

SIN(X) 



CLA 

=0202 

PACK IN EXPONENT FUR BIT 2 SCALE 


LLS 

0 

BUT BE SURE AND RETaIN SIGN OF 


LRS 

8 

COS(X)... THE RESULT IS AN 


XCA 


UNNORMALIZED FL . PT. NO. (IN ACC.) 


FAD 

= 0 

NORMALIZE COS(X) AND SAVE 


STO 

COS(X) 

THE FLOATING POINT RESULT 


LDQ 

SIN(X) 



CLA 

=0202 

PACK IN EXPONENT FOR BIT 2 SCALE AND 


LLS 

0 

MAKE SURE THE SIGN OF SIN(X) IS 


LRS 

8 

RETAINED. . .THE RESULT IS AN 


XCA 


UNNORMALIZED FL. PT. NO. (IN ACC.) 


FAD 

=0 

NORMALIZE SIN(X) AND 


XCA 


PLACE RESULT IN MQ. 


CLA 

COS(X) 

PLACE COS(X) IN ACC. 

EXIT1 

AXT 

** » 4 



AXT 

** , 1 



AXT 

**, 2 



TRA 

EJECT 

1,4 


DT 

DEC 

0 * 0245432 999B 1 

(PI/2)*(l/64)*0. 999984 (SCALED AT BIT 1) 
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Printout B-l (contd) 


DT0V2 

DEC 

0*0122716 49 9B1 

{ P I / 2 ) =M 1/64) /2*0.999984( SCALED AT 

BIT 1) 

10VDT 

DEC 

40.7436654 

1/DTHETA (FLOAT. PT. ) 


DX 

PZE 


DELTAX 


SIGN 

PZE 

** 

SIGN OF X 


T I 

PZE 

** 

COSINE OF DELTAX (SCALED AT 

BIT 1) 

T 2 

EQU 

DX 

SINE OF DELTAX (SCALED AT 

BIT 1) 

F I 

PZE 

** 

SIN(X)*T1+C0S(X)*T2 (SCALED AT 

BIT 2) 

F2 

PZE 

** 

C0S(X)*T1-SIN(X)*T2 (SCALED AT 

BIT 2) 

S IN ( X ) 

PZE 

** 

ABSOLUTE VALUE OF SIN(X) 


COS(X) 

PZE 


ABSOLUTE VALUE OF COS(X) 



EJECT 





REM 

THE FOLLOWING TABLE CONTAINS FIXED-POINT VALUES 



REM 

SCALED AT BIT 1 



SINX 

NULL 





DEC 

0*0 1 2271 5383B1 




DEC 

0# 036807 2229B1 




DEC 

0*061 32 073 63 B1 




DEC 

0.0857973123B1 




DEC 

0.1102222073B1 




DEC 

0 *134580708561 




DEC 

0.1588581433B1 




DEC 

0.1830398880B1 




DEC 

0* 2071 1 1 3762B1 




DEC 

0.2310581083B1 




DEC 

0.2548656596B1 




DEC 

0*278519689481 




DEC 

0* 3020059493B1 




DEC 

0 • 3253 102922B 1 




DEC 

0*348418680361 




DEC 

0*371317 1940B 1 




DEC 

0*3939920 40 1B1 




DEC 

0*4164295 60 1B1 




DEC 

0*438616238 5B1 




DEC 

0.4605387110B1 




DEC 

0* 48 2 1 837721 B1 




DEC 

0*50353838 37 B1 




DEC 

0.5245 89 6827B1 




DEC 

0.5453249884B1 




DEC 

0 . 56573 18 108 B 1 




DEC 

0.5857978575B1 




DEC 

0.6055110414B1 




DEC 

0.6248594881B1 




DEC 

0 *643831 5429 B 1 




DEC 

0.6624157776B1 




DEC 

0.6806009978B1 




DEC 

0.6983762494B1 




DEC 

0*7157 30 8253 B1 




DEC 

0.7326542717B1 




DEC 

0.7491363945B1 




DEC 

0*765167265661 




DEC 

0 • 7807372286B 1 




DEC 

0.7958369046B1 




DEC 

0 • 810457 1983B 1 




DEC 

0*8245893028B1 




DEC 

0.838 22 47056B1 




DEC 

0.8513551931B1 




DEC 

0.8639728561B1 




DEC 

0*876070094281 




DEC 

0.8876396204B1 




DEC 

0*898674465781 




DEC 

0 • 9091679831 B 1 




DEC 

0*91911 38517B1 




DEC 

0. 9 2 85060805 B1 




DEC 

0 *9373390 1 19B1 




DEC 

0.9456073254B1 




DEC 

0*95 33060 404B1 




DEC 

0*9604305194B1 




DEC 

0.9669764710B1 




DEC 

0.9729399522B1 




DEC 

0*97831737 07 B1 




28 


JPL TECHNICAL REPORT 32-1430 








Appendix C 

Tilted Reflector Program 


The tilted reflector program was written for an IBM 
1620 computer and outputs a deck of cards containing 
the quantities (see Eq. 12) ; 


a o (0). 

a tie), 

a 2 (d), 

ddo 


da 2 


d6 ’ 

~W 


The input data are tabular values describing a 

surface of revolution, as illustrated in Fig. C-l. Also 
specified are the translations to the new origin XF and 
ZF; a coordinate rotation is specified by the intersection 
point ZROT. 

The procedure used is as follows: 

(1) By the use of a straightforward coordinate trans- 
formation, the data 

p'(Py 0 deg) = p'(6>',90 deg) 

— p'(^,180 deg) 


are mapped into the new system for each input & 
value. Note that the cuts <£' = 0 or 180 deg map 
into <j> = 0 or 180 deg, but in general <j/ m 90 deg 
maps into a cut which will be identified by <£ 9O (0). 
The result of this step is data 

p(0fi deg), 90 ), p(0 , 180 deg) 

that are known at a set of 0 values that is in gen- 
eral different for each </> cut, and different from the 
set of desired output values. The remaining steps 
(2) through (4) are performed for each desired 
output value 0 o . 

(2) A cubic interpolation polynomial in 0 is fitted 
locally at four points to the functions p(0 , 0 deg), 
p(6,<f > 90 ), p(0,18O deg) and </> 9O (0). These polynomials 
and their derivatives are evaluated at 0 = $ 0 to 
yield <£ 9 o(0o) and dfoJdO at 6 = 0 o ; and p(^ 0 ,0 deg), 
/>( 0 o,</> 9 o), p(<9 0 ,180 deg) and dp/dO at <f> = 0, </> 90 , 
and 180 deg. 
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(3) Given the assumed form 


XF, ZF, ZROT see Fig. C-l 


p(0,0) ~ a o (0) + a x (0) cos 0 + a 2 (0)cos20 
it follows that 


>(00, 0) 


'1 

1 1 " 

do(0o) 

p(0o,4>9o) 

= 

1 

COS <f>90 cos 2(^90 

a x (0 o ) 

_p{do, 180)_ 


_1 

-1 1 

_^2(^o)_ 


T, R 0',// (Note that 6 r must be monotoni- 
cally decreasing with its index.) 

9999.0 any number > 400 (A last card indi- 
cator. NMAX must be < 200. If the 
input consists of more than 200 cards, 
the program will type out an error 
message; the program can be reset to 
read title card.) 


By the use of 0 9o (<9o) from step (2), the program 
then inverts the matrix and obtains a o (0 o )> #i(0o), 
and a 2 (0 Q ). 

For the matrix to be well-conditioned, it is suffi- 
cient that 090 ~ 90 deg. If 0 9o lies outside the range 
of 60 to 120 deg the program sets a 2 (d) — 0 and 
solves the two-dimensional matrix problem for 
a o (0) and a x (0). 


NMAX number of 0 values desired 
TOUT desired 0 values 

The program prints out the input data (p',0') and 
transformed data (p,0) at 0 = 0, 180 deg, and 0 9o . Then 
p and dp/d6 are printed out at the desired 0 values 
(Printout C-l). 

The Fourier coefficients and their derivatives are 
punched but not printed. 


(4) Also, from the assumed form, 


dp __ da 0 


do 


d$ 


+ 


da t 


dd 2 


de cos cos 2<^ 


= —a x (0) sin 0 “ 2 a 2 (0) sin 2 0 

<70 

By the use of the relationship 

dp __ dp __ dp 00 
^d0~^d0 00 ~d0 

— + [a x (e) sin 0 + 2a 2 (0) sin 20] ~ 

and the values of dp/dO and 00 9O /00 at 0 = 0 o 
obtained in step (2), values of dp/dO are obtained 
at 0 = 0 deg, 09o, and 180 deg. Then the same 
matrix inverse obtained in step (3) is used to yield 
dcio/dO, daJdB and da 2 /d$. 


The input data for the tilted reflector program are 
shown in Table C-l; the parameters are defined below: 

TITLE any alphanumeric statement 


Table C-l . Input data for the tilted reflector program 


Card 

Parameters 

Format 

1 

TITLE 



20A4 

2 

XF ZF 

ZROT 


8F10.0 

3 

TO) 

R(D 


20X / 2F10.0 


T(2) 

R(2) 



NMAX + 2 

T{NMAX) 

R(NMAX) 


20X, 2F10.0 

NMAX + 3 

9999.0 



20X, FI 0.0 

NMAX + 4 

MMAX 



1015 

NMAX + 5 

TOUT(1) TOUT(2) • 

■ TOUT{8) 

8F10.0 


TOUT{9) 

TOUT(IO) • 





• TOUT(MMAX) 

8 F 1 0.0 
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o o o n o o o o o ooo ooo 


Printout C-l . Computer printout of the tilted reflector program 


C TILTED REFLECTOR PROGRAM 

C 

DIMENSION T(200),R(200) 

DIMENSION T I T LE { 20 ) 

DIMENSION R0( 200) , To (200) , RE { 200 ) , T E ( 200 I , RN ( 200 ) , TN ( 200 ) , PN ( 200 ) 
DIMENSION TOUT ( 70 ) 

DTR=0. 017453293 
1 READ 1001, TITLE 
PUNCH 1001, TITLE 
PRINT 2001, TITLE 
READ 1002,XF,ZF,ZR0T 
PRINT 2002, XF,ZF,Z ROT 
DZ=ZROT-ZF 

CALL ATANXY(DZ,XF,BETA) 

BEDE G= BETA/ DTR 
PRINT 2003 , BEDEG 
SlB=SINF(BETA) 

COB=COSF ( BETA ) 

READ IN SHAPED SUBREFLECTOR 

N = 0 

10 N=N+1 

READ 1003 , T ( N ) , R ( N ) 

T(N )=T (N)*DTR 
IF(T(N>-4. 0)15,20, 20 
15 IF(N-200) 10,99,99 
20 NMAX=N- 1 

PRINT 2004, NMAX 

TRANSFORM DATA TO NEW SYSTEM 

CALL TRANQ(R,T,RO,TO,PN,0„0, 1 . 0 , S I B ,C OB , XF , ZF , NMAX ) 

DO 24 N=1 , NMAX 

24 T0(N)~T0(N)#C0SF(PN(N) ) 

CALL TRANQ(R,T,RE,TE,PN,0.0,-1.0,SIB,COB,XF,ZF,NMAX) 

DO 25 N=1 , NMAX 

25 TE(N)=-TE(N)^COSFtPN(N) ) 

CALL TRANQ(R,T,RNtTN,PN,1.0, 0. 0, SIB, COB, XF,ZF, NMAX) 

PRINT 2008 
DO 26 N=1 , NMAX 
TOO=TO(N) /DTR 
TNN=TN(N)/DTR 
TEE=TE(N)/DTR 
TTT=T(N)/DTR 
PNN=PN(N)/DTR 

26 PRINT 2009, TTT, R(N) »TOO,RO(N) , TEE, RE(N) ,TNN,RN(N) ,PNN 

READ DESIRED THETA VALUES 

READ 1004, MMAX 
I F ( MMAX- 75)29,29,99 
29 READ 1002, ( TOUT ( N ) , N= 1 , MMAX) 

NLAS TG=NM AX-2 
NLASTE=NMAX-2 
NLASTN=NMAX-2 
PRINT 2007 
DO 30 MM= 1 , MMAX 
TTO=TOUT(MM)*DTR 

INTERPOLATE TO DESIRED THETA VALUES 

CALL QBURP < TO ,R0 » TTO , RRO ,DRO,NLASTO, NMAX ) 

CALL QBURP (TE, RE, TTO, RRE, DRE, NLASTE , NMAX ) 

CALL QBURP (TN,RN, TTO, RRN, DRN, NLASTN, NMAX) 

CALL QBURP ( TN,PN, TTO, PPN,DPN, NLASTN, NMAX ) 

PNN=PPN/DTR 

FIND FOURIER COEI FF IC I ENTS 
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Printout C-l (contd) 


CALL FKQSF(RRO,RRE,RRN,PPN,AO,A1,A2) 

DRN=DRN+( A I* S IMF ( PPn ) + 2 . 0*A 2*S I NF < 2 . o*PPN ) ) *DPN 
PRINT 200 5, TOUT (MM) , RRO , DRO T RRE , DRE , RRN , DRN , PNN , DPN 
CALL FKQSF < DRO » DRE , DRN , PpN , DAO , DAI , DA2 ) 

PUNCH 2006, TOUT (MM) ,A0, Al, A2 , DAO , DA 1 , DA2 
30 CONTINUE 
GO TO 1 
C 

99 TYPE 2010 
PAUSE 
GO TO 1 
C 

1001 FORMAT (20A4) 

1002 FORMAT ( 8F10.0 ) 

1003 FORMAT (20X,2F10.0) 

1004 FORMAT ( 1015 ) 

2001 FORMAT! 1H1,20A4) 

2002 FORMAT ( 28H COORDINATES OF NEW ORIGIN/ 

*7H X= , E 15 • 8/ 

*7H Z=,E15.8/ 

*2 6H POINT OF ROTATION ZROT= , El 5. 8 ) 

2003 FORMAT ( 23H ROTATION ANGLE BETA=,F8.4,9H DEGREES) 

2004 FORMAT ( 31 H NUMBER OF INPUT POINTS NMAX=,I4) 

2005 FORMAT ( F10.5,3(F14.5»F12.5) ,2F11.5) 

2006 FORMAT ( 7F 10 • 4 ) 

2007 FORMAT! 1H1 ,47H R AND ITS PARTIAL DERIVATIVE WITH RESPECT TO, 

*3iH theta at desired output points// 
*20X5HPHI=0,21X7HPHI=180,21X16HPHI AS TABULATED/ 

*10H THETA ,3<8X,1HR,10X,7HDR/DT ) , 4X3HPHI , 7X7HDPHI /DT ) 

2008 FORMAT ( //52H SURFACE DATA IN INPUT AND TRANSFORMED COORDINATES// 
*7X , 10H INPUT DATA, 11X»15HPHI= 0 DEGREES , 10X , 15HPH I =1 80 DEGREES, 15X 
*,16HPHI AS TABULATED/ 

*4(8H THETA, 8X, 1HR,8X) ,4H PHI) 

2009 FORMAT! F 10 . 5 , F 1 2 . 5 , 3 ( F 1 3. 5 , FI 2 . 5 ) , FI 1.5) 

2010 FORMAT (41H INPUT DATA EXCEEDS DIMENSION STATEMENT/ 

*42 H RELOAD DATA AND PUSH START TO TRY AGAIN) 

END 

SUBROUTINE QBURP ( T , R , TT , RR , DR , NLAST , NMAX ) 

DIMENSION T( 2) ,R(2) 

DIMENSION 8(4,5) , A ( 4 ) 

THIS SUBROUTINE NUMERICALLY INTERPOLATES AND DIFFERENTIATES THE 
FUNCTION R ( T ) BY FITTING A CUBiC TO R ( T ) AT FOUR POINTS 

FIND NLAST SUCH THAT T(NLAST-l) AND T ( NLAST ) STRADDLE THE 
DESIRED VALUE TT. IT IS ASSUMED THAT THE TABULAR VaLUES OF T 
DECREASE WITH INCREASING N 
9 I F (T ( NLAST ) -TT ) 10,10,20 
10 IF(NLAST-2)20,20,15 
15 NLAST =NLAST~1 
GO TO 9 
20 CONTINUE 

SOLVE JHE SYSTEM OF FOUR EQUATIONS TO DETERMINE THE COE I FF I Cl ENTS 

ESTABLISH MATRIX 
DO 30 1=1,4 
NK=NLAS T+I “2 
B ( I , 5 ) =R ( NK ) 

DO 30 J= 1 , 4 
JJ= J~1 

30 B(I»J)=T(NK)**JJ 

DIAGONALIZE SYSTEM 
DO 50 K= 1 , 3 
KK=K+ 1 

DO 50 I=KK, 4 
D=B ( I , K ) /B ( K, K ) 

DO 50 L=K , 5 

50 B ( I , L ) = B ( I , L ) -0*B ( K, L ) 
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Printout C-l (contd) 


c 


SOLVE BY BACK SUBSTITUTION 


A(4)=B(4,5)/B(4,4> 
DO 70 11=1*3 


1=4-11 


L=5-II 


DO 60 J = L » 4 

60 B(I ,5)=B( I,5)-B(I t J ) *A ( J) 
70 Am=B(I,5)/B(I,I) 


EVALUATE CUBIC And its derivative at interpolation POINT Tt 


RR=A ( 1) +A( 2)*TT+A(3) *TT*TT+A ( 4 ) *TT*TT*TT 

DR=A(2)+2.0*A(3)*TT+3.0*A(4)*TT*TT 

RETURN 

END 

SUBROUTINE FKQS F ( 0 , E * FN , P , AO , A1 , A2 ) 


this Subroutine computes the cgeifficients for the truncated 

FOURIER SERIES R( THETA ) »A0+A1*C0S < THETA )+A2*C0S ( 2*THETA ) 

0=R(ZER0 RADIANS) 

E=R (PI RADIANS) 

FN=R ( P RADIANS) WHERE P SHOULD LIE BETWEEN 1.05 AND 2.09 

IF(P-1. 05)100, 20, 20 
20 I F ( 2.0 9- P )100,200, 200 
100 A0=(0+E)/2.0 
A 1= ( G-E ) / 2. 0 
A2=0. 

RETURN 

200 C = COS F ( P ) 

C2= COSF ( 2 • 0#P ) 

D=2.0*(C2“1.0) 

A0= (C+C2) *Q-2.0*FN+(C2-C )*E 
A0=A0/D 

Ai= (C2-1. 0)*0+ (-C2+ 1.0) *E 
A1=A1/D 

A2= (-C-1«0)*0+2.0*FN+(C-1.0)*E 

A2-A2/D 

RETURN 

END 


SUBROUTINE TR ANQ ( R , T , RR , TT , PP , S I ,C0 , S I B ,COB , XF , ZF , NMAX ) 

THIS SUBROUTINE PERFORMS A COORDINATE TRANSLATION AND ROTATION 

DIMENSION R(2),T(2),RR(2) , TT ( 2 ) » PP ( 2 ) 

C 

C CONVERT TO RECTANGULAR COORDINATES 

DO 10 N=1,NMAX 
X=R(N)*SINF(T(N) )*C0 
Y=R (N)*SINF(T<N) ) -SI 
Z=R(N)*cOSF(T(N) ) 

C 

C PERFORM TRANSLATION 

X=X -XF 
Z=Z-ZF 
C 

C PERFORM ROTATION 

D=X 

X=X*COB+Z*SIB 
Z=— D*S I B+Z^COB 
C 

C CONVERT TO POLAR COORDINATES IN NEW SYSTEM 

RR(N)=SQRTF(X*X+Y*Y+Z*Z ) 

D=SQRTF(X*X+Y*Y) 

CALL AT ANXY ( Z , D , TT ( N ) ) 
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Printout C-l (contd) 


10 CALL ATANXY { X » Y * PP ( N) ) 
RETURN 
END 

SUBROUTINE ATANXY(X,Y,PHI ) 
P I =3 • 1415927 
I F ( X ) 10,20,30 
10 PHI=ATANF(Y/X)+PI 
GO TO 40 

20 IF(Y)21 ,22,23 

21 PHI=-PI/2.0 
GO TO 50 

22 PH I =0 • 0 
GO TO 50 

23 PHI=PI/2.0 
GO TO 50 

30 PHI=AtANF (Y/X) 

GO TO 50 

40 I F ( PHI -P I ) 50 » 50 » 41 

41 PHI=PHI-2.0*PI 
50 RETURN 

END 
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Appendix D 

Scattering and Spherical-Wave Expansion Programs 


The scattering program consists of the following pro- 
grams* (see Fig. D~l): 

MAIN asymmetrical scattering program 

FIELDS subroutine that provides H* on surface 
(see Section IV) 

SURF subroutine that provides p and its deriva- 
tives (see Section III) 

SETUP subroutine that establishes integration and 
output grids, and tables of precomputed 
trigonometric functions 

PATHL subroutine that established the path length 
term, Eq. (17) 

FINT numerical integration subroutine (see Sec- 
tion V) 

LEGEND Legendre functions subroutine 
SPHANK spherical Hankel function subroutine 

EXPJX fast trigonometric subroutine (see Appen- 
dix B) 

VECTOR subroutine that converts from rectangular 
to polar coordinates 

ADJUST subroutine that normalizes angles to the 
range — 180 to 180 deg 

PRTIM subroutine that prints out execution times 

The spherical-wave expansion program (orthogonality 
method) consists of the following subprograms (see 
Fig. D-2): 

MAIN spherical-wave expansion program (or- 
thogonality method) 

MULT a general purpose matrix multiplication 
subroutine 

LEGEND Legendre function subroutine 

VECTOR subroutine that converts from rectangular 
to polar coordinates 

PRTIM subroutine that prints out execution times 

*A11 the programs presented in this section, except SOLVE and 
PRTIM, were developed and programmed by the author. The 
matrix inversion package SOLVE was developed by Richard 
Hanson of the Computation and Analysis Section of JPL; PRTIM, 
which prints out the actual computation time during program 
execution, was developed by Charles Lawson, also of the Com- 
putation and Analysis Section. 


The spherical-wave expansion program (linear equa- 
tion method) consists of the following subprograms (see 
Fig. D-2): 

MAIN spherical-wave expansion program (linear 
equation method) 

SOLVE a package of subroutines to solve a set of 
linear equations 

LEGEND Legendre function subroutine 

PRTIM subroutine that prints out execution times 

The input data required by the main scattering pro- 
gram are shown in Table D-l; the parameters are as 
follows: 

TITLE any alphanumeric statement 

PC propagation constant, %r/\ 

XT, YT, ZT translations to reference point for 
phase pattern of output data (nor- 
mally the expected phase center of the 
scattered pattern) 

SCALE a scale factor for the output fields; if 
this is left blank or if input is zero, 
program sets SCALE = 1.0 

JMAX number of 0 values desired, JMAX 
<181 

TT1 initial 0 value 

DTT increment for 0 values, DTT > 0 

TT(1), TT(2), • • • explicit list of 0 values 

KMAX number of d> values desired in output, 
KMAX < 8 

PP1 initial $ value 

DPP increment for <£ values, DPP > 0 

PP(1), PP(2), • * • explicit list of <£ values 

NCI number of integration grids, NCI < 5 

MM (I) number of 0 values in the Ith integra- 
tion grid, MM < 15 

T1 initial 0 value 

DT increment for 0 values, DT > 0 
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(a) ORTHOGONALITY VERSION 



(b) LINEAR EQUATION VERSION 


^ START ^ 





READ INPUT PATTERN 





INVERT SYSTEM OF LINEAR 
EQUATIONS TO OBTAIN WAVE 
COEFFICIENTS; OUTPUT 
RESULTS 



Fig. D-2. Block diagram of spherical-wave 
expansion program 


T(I,1), T(I,2), ■ • • 
NN(I) 


PI 

DP 


explicit list of 0 values for Ith grid 

number of $ values in the Ith integra- 
tion grid, NN < 91 

initial <f> value 

increment for <j> values 


P(I>1)> P(I>2), 

(LAST) 


explicit list of < p values for Ith grid 

a last case indicator; punched as 
shown in columns 1 through 6 


1 THIS IS A JMAX X KMAX LOOP 

2 THIS STEP INTERNALLY CONTAINS 

SEVERAL MMAX x NMAX LOOPS 

Fig. D-l . Block diagram of scattering program 


In all cases, a set of angle values may be specified by 
an initial value and a positive increment, or by an ex- 
plicit list of values (not necessarily equally spaced). All 
angles are in degrees. If a positive increment is given, the 
data cards containing explicit values must be omitted. If 
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Table D-L Input for MAIN scattering program 


Table D-2. Input for spherical-wave FIELDS subroutine 


Parameters 


TITLE 

1 3A6 

PC XT YT ZT SCALE 

5F10.0 

JMAX TT1 DTT 

15,2 F 1 0.0 

TT{1) TT (2) • • • 

8F10.0 

• ■ • TT(JMAX) 

8F10.0 

KMAX PP1 DPP 

15,2F1 0.0 

PP{1) PP(2) ♦ • * 

8F10.0 

• ■ ■ PP(KMAX) 

8F10.0 


IJC 

INV7 I 

13 

MM(1) T1 DT 

15,2 F 1 0.0 

T(l,1) T(l,2) • ■ • 

8F10.0 

• • • T{1 ,MM(1 )) 

8F10.0 

NN(1) PI DP 

I5,2F 1 0.0 

P(1,1) P{1,2) • • * 

8F10.0 

• • • P(1,NN(1)) 

8F10.0 


MM(NGI) T1 DT 
T(NG1 ,1 ) T(NG1 ,2) • • • 

• • * T(NG1,MM(NGT)) 
HN(NGl) PI DP 
P(NG1,1) P(NG1 ,2) • ■ • 

* ♦ • P(NG1 ,NN(NG1)) 
[Input for FIELDS Subroutine] 
[Input for SURF Subroutine] 
TITLE 

PC XT YT ZT SCALE 


(Data for further cases, 
same as above) 


I5,2F1 0.0 
8F10.0 
8F10.0 
I5,2F1 0.0 
8F1 0.0 
8F10.0 


the data cards containing explicit values are included, 
the increment must be written in as identically zero. 

The input data required by the spherical-wave version 
FIELDS subroutine are shown in Table D-2; the param- 
eters are as follows: 

TITLE any alphanumeric statement 

LMAX maximum mode order, < 60 

MCOMP order of azimuthal variation 


Card 

Parameters 

Format 

1 

TITLE 

1 3A6 

2 

TITLE 

1 3A6 

3 

LMAX MCOMP 

215 

4 

1 A(1 ,1) A(1 ,2) B(1 ,1) B{1 ,2) 

15,2 El 7.8, 
2X,2E1 7.8 

LMAX + 2 

LMAX A(LMAX,1) A{LMAX,2) 
B(LMAX,1 ) B(LMAX,2) 

• 


A(N,1), A(N,2) real and imaginary parts of TE MCOMP, N 
wave coefficient (see Appendix A and 
following program) 

B(N,1), B(N,2) real and imaginary parts of TM MCOMP, N 
wave coefficient (see Appendix A and 
following program) 

The input data for the far-field version FIELDS sub- 
routine are shown in Table D-3; the parameters are as 
follows: 

TITLE any alphanumeric statement 
JMAX not used 
JO not used 

JIN number of input points, JIN < 181 
IC1 < 0 for input in dB; > 0 for input in V 
IC2 > 0 to set phase pattern zero everywhere 
MCOMP order of azimuthal variation 
PSI polar angle 0 (see Fig. 1) 

E I E e (<9, 90 deg) I, in V or dB 


EP argument {E e ((9, 90 deg)}, in deg ) S ee 


H | E<p(6 , 0 deg) |, in V or dB 
HP argument {£0(0,0 deg)}, in deg 


Eq. (14) 


The input data for the SURF subroutine are shown in 
Table D-4; the parameters are as follows: 

TITLE any alphanumeric statement 

TIN input 6 value (see note below) 
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Table D»3. Input for far-field FIELDS subroutine 


Card 

Parameters 

Format 

1 

TITLE 

72H 

2 

JMAX JO JIN IC1 IC2 
MCOMP 

615 

3 

PSI(l) £(1 ) EP(1 ) H(1 ) HP(1) 

5F10.6 

JIN + 2 

PS! (JIN) E(JIN) EP(JIN) H(J1N) 
HP(JIN) 

5F10.6 


AO, Al, A2 a 0> a ly a 2 (see Eqs. 12 and 13, and 
Appendix C) 

DAO, DAI, DA2 (see Eqs. 12 and 13, 

and Appendix C) 

The TIN values must agree with the specified inte- 
gration grid values. (If they disagree, an error message is 
printed but the program continues to run.) The values 
for each integration grid are stacked in order. 

Computer time is nearly proportional to the quantity 
P - £ NMAX X MMAX 

NGl integration 
grids 

For an IBM 7094, the time can vary by a factor of 2 
depending on which model is being used; time for a 
slower model is given by 

t = P X [JMAX X KMAX] 8.4 X Kh 5 min 

plus 

t = PX [LMAX] 2 5 X 10- 6 min 
if the spherical-wave version FIELDS subroutine is used. 

The input data required by the spherical-wave expansion 
program (orthogonality version) are shown in Table D-5. 
All of these parameters have been previously defined, 
except for JMAXO = 180/A# + 1, where A 0 is the 
desired output increment, and JOUT is the number of 
output values starting with 9 = 0. The program evaluates 
the spherical-wave expansion at p = oo and prints and 


Table D-4. Input for SURF subroutine 


Card 

Parameters 

Format 

1 

TITLE 

12A6 

2 

TIN(l) A0(1 ) Al (1 ) A2(l ) DA0(1) 
DAI (1 ) DA2(1 ) 

TIN(2) A0{2) Al (2) A2(2) DA0(2) 
DAI (2) DA2{2) 

7F10.4 


TIN(MM) AO(MM) Al (MM) A2(MM) 
DAO(MM) DAI (MM) DA2(MM) 
TlN.fi ) A0(1) Al (1 ) • ■ • 

7F10.4 


punches the results at a set of output values determined 
by these parameters. 

The input data required by the spherical-wave expan- 
sion program (linear equation version) are identical except 
that the first two cards and the last card are omitted. 

Both of these programs use less than 1 min of computer 
time to evaluate 60 TE and 60 TM wave coefficients. 

Two computer printouts are included in this appendix: 
the first (Printout D-l) is a computer printout of the asym- 
metrical scattering program; the second (Printout D-2) 
is a printout of the spherical-wave expansion program. 


Table D-5. Input for spherical-wave expansion program 


Card 

Parameters 

Format 

1 

TITLE 

13A6 

2 

MCOMP LMAX 

215 

3 

TITLE 


4 

JMAX JO JIN IC1 IC2 
MCOMP 

615 

5 

PSI(l) E(1 ) EP(1 ) H{1 ) HP(1 ) 

5F10.0 

JIN 4- 4 

PSl(JIN) E(JIN) EP(JIN) H(JIN) 
HP(JIN) 

5F10.0 

JIN 4* 5 

JMAXO JOUT 

215 
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Printout D-1. Computer printout of asymmetrical scattering program 

(a) Main program 

$ I BFTC MAIN 

c Asymmetrical scattering program 

C PROGRAM COMPUTES SCATTERING FOR PERFECTLY CONDUCTING SMOOTH 

C SURFACE OF ARBITRARY SHAPE A . LUDWIG 1-15-69 

C 
C 
C 

COMMON/GRID 1/SI T( 15) ,COT( 15) ,$iP{ 181) ,COP( 18 1 ) , T ( 15 ) , P (18 1 ) 
COMMON/GRI D2/SITT(181),CQTT(181),SIPP( 8) ,COPP ( 8 ) , TT ( 18 1 ) , PP { 8 ) 

DIMENSION TITLE! 13 ) 

DIMENSION F { 15, 91),FT( 15, 9l>tFP( 15, 91),GAM( 15, 91) 

INTEGER TITLE 

COMPLEX HR ( 15, 91),HT( 15, 91),HP( 15, 91) 

COMPLEX ETT (181, 8),EPP(181, 8) 

COMPLEX ETTO , EPPO 
COMPLEX A ( 15 , 91,3) 

COMPLEX STOT ( 3 ) 

COMPLEX T1,T2,T3 
COMPLEX TET , TEP 

EQUIVALENCE ( A ( 1 , 1 , 1 ) , HR ) , ( A { 1 , 1 , 2 ) , HT ) , ( A ( 1 , 1 , 3 ) , HP ) 

EQUIVALENCE (GAM, FT) 

DATA LC/ 6H( LAST ) / 

C START TIMING ROUTINE 

CALL PRTIM1 
READ(5, 100DTITLE 
1 WRITE(6, 2001) TITLE 

READ(5» 100 2) PC, XT ,YT ,ZT, SCALE 
IF ( SCALE) 22,21,22 

21 SCALE=1.0 

22 CONTINUE 
WRITE(6,2003)PC 

READ IN GRID DATA AND ESTABLISH OUTPUT GRID AND FIRST INTEGRATION 
GRID 

CALL SETUP (NG1, 1 , JMAX , KMAX » MMAX , NMAX ) 

CALL PRTIM2 
WRI TE ( 6, 2002 ) TI TLE 

READ IN EXCITATION FIELD DATA 
CALL FIELDS 
CALL PRTIM2 

BEGIN LOOP FOR INTEGRATION GRID SEGMENTS 
DO 90 J = 1 » JMAX 
DO 90 K=1,KMAX 
ETT(J,K)=<0. 0,0.0) 

90 EPP(J,K)=(0. 0,0.0) 

100 WR ITE (6 , 2002 ) TI TLE 

ESTABLISH SURFACE PARAMETERS ON INTEGRATION GRID 
CALL SURF(MMAX,NMAX,F,FT,FP) 

DO 410 M=1,MmAX 
DO 410 N=1,NMAX 
410 F(M,N)=PC*F(MtN) 

Establish excitation fields on integration grid 
call fieldkmmax,nmax,hr,ht,hp»f) 

combine surface and field data to determine complex vector a 

DO 400 M=1,MMAX 
DO 400 N=1,NMAX 

T1 = FT (M »N)#SIT(M)*HP(M,N)-FP(M,N)*HT{M,N) 
T2-FP(M,N)*HR(M,N)+F(M,N)*SIT(M)*HP(M,N)/PC 
T3=~FT{M,N )*S1T(M)*HR(M,N)-F(M,N)*SIT(M)*HT (M,N)/PC 
A ( M , N , 1)=T1*SIT(M)*C0P(N)+T2*C0T(M)#C0P(N)-T3*SIP(N) 

A(M,N,2 )=T1*$IT(M)*SIP(N)+T2*C0T(M)*SIP(N)+T3*C0P(N) 
A(M,N,3)=T1*C0T( M)-T2*SIT(M) 

400 CONTINUE 

WRITE(6, 2007)1 
CALL PRTIM2 

BEGIN OUTPUT GRID LOOP 
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Printout D-l (contd) 


WRITE(6, 2006)1 
C 

DO 500 K=1»KMAX 
TO-PP(K) /0«017453293 
WRITE(6,2011)T0 
DO 500 J= 1 , JMAX 
C 

c Establish phase/path length parameter on integration grid 

Call PATHL(F,JtKtMMAX,NMAX»GAM) 

c 

C PERFORM INTEGRATION 

CALL F INT(T,P,A,GAM,MMAX-1,NMAX~1,ST0T) 

C 

C CONVERT RESULT TO PoLAR COORDINATES AND SUPERINPOSE FIELDS 

C 

ETTO=COTT! J)*(STOT! 1 ) *COPP t K ) +STOT ( 2 ) *S I PP ( K ) ) -STOT < 3 ) *SI U ( J) 
EPPO-STOTI 2)^C0PP(K)-ST0T{ 1)*SIPP(K) 

ETT0=-<0.0, 1.0 )*PC/ 6. 283 1854* ETTO 
EPPO=-{ 0.0 , 1.0) *PC/6. 2831 854*EPP0 
TO=TT< J)/0. 017453293 
A1=REAL !ETTO) 

A2= A IM AG { ETTO ) 

A3=REAL ( EPPO ) 

A4=AIMAG! EPPO} 

CALL VECTOR! A1 , A2 , ETAMP , ETPHI ) 

CALL VECTOR! A3, A4 , E PAMP , EPPH I ) 

WRITE! 6, 2012) TO, ETAMP, ETPHI ,EPAMP,EPPHI 
ETT (J,K)=ETT!J»K)+ETTO 
EPP(J ,K )=EPP(J,K)+EPPQ 
500 CONTINUE 

CALL PRTIM2 
C 
C 

C IF MORE INTEGRATION GRIDS REMAIN LOOP BACK 

C 

I F ( I -NG 1)600,700,700 
600 CALL RESET! I ,MMAxtNMA X ) 

GO TO 100 
700 1=1+1 

C ESTABLISH DIRECT RADIATION ON OUTPUT GRID 

WRI TE ( 6, 2009 ) 

DO 750 K=1,KMAX 

T0=PP(K)/0*O17453293 

WRITE!6,2011)T0 

DO 750 J= 1 , JMAX 

CALL FIElD2( J,K,ETTO,EPPO) 

TO=TT(0)/0. 017453293 
A1=REAHETT0) 

A2=AIMAG! ETTO) 

A3=RE AL ! EPPO ) 

A4= AIMAG! EPPO ) 

CALL VEcTORfAl, A2 , ETAMP, ETPhI ) 

CALL VECT0R(A3,A4,EPAMp,EPPHl) 

WRITE (6, 201 2) TO, ETAMP, ETPHI ,EPAMP, EPPH I 
ETT! J,K)=ETT! J,K)+ETTO 
EPP(J,K)=EPP!J,K)+EPPO 
750 CONTINUE 

CALL PRTIM2 
C 

C TRANSLATE PHASE CENTER, SCALE AMPLITUDES, AND OUTPUT TOTAL FIELDS 

C 

PUNCH 1001, TITLE 
WRITE!6,2002)TITLE 
WRITE!6,2010)XT,YT,ZT,SCALE 
DO 760 K= 1 , KMAX 
T0=PP(K)/O. 017453293 
WRITE!6,2011)T0 
DO 760 J=1 , JMAX 
T0=TT(J)/O. 017453293 
A 1 =RE AL ! ETT ! J , K ) ) 
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Printout D-l (contd) 


A2=AIMAG(ETT(J,K) ) 

A3=REAL(EPP( J,K) ) 

A4=AIMAG(EPP( J,K) ) 

CALL VECTOR (AItA2* ETAMp, ETPH I ) 

CALL VECT0R(A3,A4,EPAMP,EPPHl ) 

DP=XT*SITT( J)*CQPP(K)+YT*SITT( J)*SIPP(K)+ZT*(C0TT(J)+1.0) 

DP=DP*PC*57. 29578 

ETPHI=ETPHI“DP 

epphi=epphi-dp 

CALL ADJUST { ETPHI ) 

CALL ADJUST(EPPHI ) 

ETAMP=ETAMP*SCALE 

EPAMP=EPAMP*SCALE 

WRITE < 6,2012) TO, ETAMP, ETPHI ,EPAMP,EPPHI 
PUNCH 2013, TO, ETAMP, ETPHI ,EPAMP,EPPHI 
760 CONTINUE 

CALL PRTIM2 


800 


1001 

1002 

2001 

2002 

2003 

2004 
2006 
2007 

2009 

2010 


TEST FOR ANOTHER DATA 
READ(5»1001)TITLE 
IF (TITLE! 1)-LC) 1,800,1 
CONTINUE 
WRITE (6,2004) 

STOP 

FORMAT ( 13A6 ) 

FORMAT (10F10.0) 

FORMAT ( 32H 1 
FORMAT (5H1 
FORMAT ( 26H 
FORMAT ( 18 HI 
FORMAT ( 29H0 
FORMAT ( 30H1 
FORMAT (40H1 
FORMAT ( 63H 
1 FIELDS/32H 


BLOCK 


LUDWIG SCATTERING PROGRAM/ /5X , 1 3A6// ) 

, 1 3A6// ) 

PROPAGATION CONSTANT^ , El 4 # 8 // ) 

END OF LAST CASE) 

SCATTERED FIELDS FROM GRID, 12) 

BEGIN INTEGRATION OVER GRID, 12) 

DIRECT RADIATION FROM INCIDENT FIELDS) 
SUPERPOSITION OF ALL GRID SCATTERED FIELDS AND 
PHASE CENTER TRANSLATED BY X=,F10,4,5H Y = ,F 


DIRECT 
10 • 4 , 


25H Z= » F 1 0 . 4/ 

33 9h amplitud e values scaled by factor of,ei5.8) 

2011 FORMAT ( 7H0 PHI=,F7.2/ 

117X,7HE THETA, 15X,5HE PHI/ 

29H THETA, 2(20H VOLTS PHASE )) 

2012 FORMAT (F9.2.FU.6»F8.2|F12.6,F8.2) 

2013 FORMAT ( F 10. 2, F 10 ,6 » F 10 • 2, F 10. 6, Flo* 2) 

END 


(b) Subroutine FIELDS (spherical-wave expansion) 


$ I BFTC FIELD 

SUBROUTINE FIELDS 
C 

C THIS SUBROUTINE EVALUATES A SPHERICAL WAVE EXPANSION TO OBTAIN 

C THE MAGNETIC FIELDS ON A SURFACE S. WAVE COE I FF IC I EnT$ ARE INPUT 

C ON CARDS. A. LUDWIG 1-14-69 

C 

COMMON/ GR I D1/SIT( 15) ,COT( 15) ,SIP( 181 ) ,COP( 181 ) , T( i5 ) , P( 181 ) 
C0MM0N/GRID2/SITT(181),C0TT(181),SIPP(8),C0PP(8),TT(181),PP(8) 
DIMENSION F ( 10 1 ) , G ( 100 ) 

DIMENSION NAME ( 13 ) 

DIMENSION A( 100,2) ,B< 100,2) 

C 

C READ IN WAVE COE I FF IC I ENTS 

C 

RE AD ( 5 , 100 1 ) NAME 
WRITE(6,2001)NAME 
READ ( 5 , 100DNAME 
RE AD (5, 1002) LMAX,MCOMP 

READ (5, 1003) (J, A (J, 1) , A ( J , 2 ) , B ( J , 1 ) , B ( J , 2 ) , L= 1 , LMAX ) 

FMC=MCOMP 

WRITE (6, 2003) MCOMP 
WRI TE(6»2002)NAME 

WRITE (6, 2004) ( J, A( J, 1 ),A(J,2),B(J,1),B(J,2) , J=1,LMAX) 

RETURN 
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Printout D-l (contd) 


c 

C ENTRY POINT FOR MAGNETIC FIELDS AT FINITE R 

C 

ENTRY FIELDHMMAX,NMAX,HR,HT,HP,R) 

DIMENSION R ( 15, 91) 

COMPLEX HR ( 15, 91),HT( 15, 9l),HP( 15, 91) 

I ENT = 1 
M~0 

1 M = M + 1 
SN=SIT (M) 

Z-COT(M) 

TOUT=T ( M ) 

GO TO 99 
C 

C ENTRY POINT FOR ELECTRIC FIELDS AT INFINITE R 

C 

ENTRY FIELd2( JO , KO , ETTO , EPPO ) 

COMPLEX ETTO, EPPO 

IENT=2 

SN=SITt(JO) 

Z^COTTUO) 

TOUT=Ty { JO ) 

c 

C ESTABLISH LEGENDRE FUNCTION TABLES FOR M-TH THETA VaLUE 

C 

99 I F ( AB S ( SN ) -*0000 1)200, 100, 100 
lOO DO 105 N=1 ,mC0MP 
105 F ( N ) = 0 

NC=LMAX+1 

CALL LEGEND(NC,MCOMP,Z,F) 

DO 110 N= 1 , LMAX 

T1=N-MC0MP+1 

T2=N+1 

G(N)=T1*F(N+1)-T2*Z*F( N) 

110 G(N)=G(N)/SN 

DO 115 N=1,LMAX 
115 F(N)=F(N)/SN 

C F ( L ) IS MULTIPLIED BY FMC LATER 

GO TO ( 300 , 500 ) , I ENT 

C SPECIAL EQUATIONS FOR TH=0 AND TH=180 DEG 

200 I F (MCOMP-1 ) 210,220,210 
210 DO 215 N=1 , LMAX 
F (N )=0 
215 G(N)=0 

GO TO ( 300 , 500 ) , I ENT 
220 DO 225 N=1,LMAX 
FN=N*(N+1) 

F(N)=FN/2.0 
225 G ( N ) =FN/2 • 0 

IF (TOUT-1. 57) 250,250, 230 
230 DO 235 N=1,LMAX,2 
F(N+1)«-F(N+1) 

235 G(N)--G(N) 

250 GO TO (300,500) ,IENT 
C 

C FOR EACH PHI VALUE THE HANKEL FUNCTIONS ARE EVALUATED AT THE 

C CORRESPONDING VALUE FOR R AND THE WAVE EXPANSION IS SUMMED 

C 

300 DO 450 N=1,NMAX 
HRR =0 
HR I =0 
H T R = 0 
HT I =0 
HP R = 0 
HP 1=0 

CALL SPHANK( 1,R(M,N) , SOR , SO I ) 

DO 400 L=1 , LMAX 
FL= L 

FN=(FL+1.0)/R(M,N) 

NC=L+1 

CALL SPHANK ( NC ,R (M,N ) ,S1R»S1I ) 
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Printout D-l (contd) 


F1R=-A { L * 1 ) #S1R+A{ L, 2 )*$1 I 
F1I=-A(L,1)*S1I-A(L,2)*$1R 
F 2 R ~ F N * ( A ( L , 1 ) *S0I + A ( L, 2 ) *SOR ) 

F2I=FN*( A(L,2)*S0I-A( L, 1)*S0R) 

F 3R = -B ( L » 1)*S0R+B( L,2)*S0I 
F3I=-BU,1)*S0I-B(L, 2)*S0R 
HRR«HRR+F2R*F(L)*Ft 
HRI=HRI+F2I*F{L)*FL 
F ( L ) ~F ( L ) *FMC 

HTR=HTR+F1R#G(L)+F2R*G(L)+F3R*F(L) 

HTI=HTI+F1I*G(L>+F2I*G(L)+F3I*F( L) 

HPR = HPR-FXR*F< L)-f2R*F(L)-F3r#G(L) 
HPl=HPI-FlI*F(U-F2I*FtU-F3I*G(L) 

S0R=S1R 
S0I=S1I 
400 CONTINUE 

HRR=HRR*SN*C0S ( FMC*P ( N ) ) 

HRI«HRI*SN*COS(FMC*P(N> ) 

HR(M,N)=CMPLX(HRR,HRI ) 

HTR=HTR*C0S( FMC*P( N) } 

HTI=HTI*cOS(FMC*P(N) ) 

HT<M,N)=CMPLX(HTR,HTI) 

HPR=HPR$$I N t FmC^P ( N ) ) 

HPI=HPI*SIN(FMC*P(N) ) 

HP(M f N)=CMPLX(HPRtHP I ) 

450 CONTINUE 

IF(M-MMAX) 1,460,460 
460 RETURN 
500 CONTINUE 
ETR=0.0 
ET I =0 
EPR = 0 
EP I =0 

DO 550 L=1»LMAX 
F(L)*F(L)*FMC 

ETR=ETR+A(L, 1)*F(L}+B(L, 1)*G( L) 

ETI=ETI+A(L,2)*FU)+B( L,2)*G(L> 

EPR=EPR+A ( L , U*G(L)+B(L, 1)*F(L) 

EPI=EPI+A< L» 2)*G(L>+B( Lt2)*F(L) 

550 CONTINUE 

ETR=ETR*SIN<FMC*PP(KO) ) 

ETI~ETI*$IN(FMC*PP(KO) ) 

ETTO=CMPLX!ETR,ETI ) 

EPR = EPR*COS ( FMC^PP < KO ) ) 

EPI-EPI*COS(FMC*PP(KO) ) 

EPPO=CHPLX (EPR,ePI) 

RETURN 

1001 FORMAT ( 1 3A6 ) 

1002 FORMAT 1515) 

1003 FORMAT( I5,2E17.8,2X,2E17.8) 

2001 FORMAT ( 1H »62H FIELD DATA INPUT IN THE FORM OF SPHERICAL WAVE CO 
1EFFICIENTS/5X,13A6> 

2002 FORMAT (5X, 13A6///20X,4HA{N) ,32X,4HB ( N ) / 

15H N,7X4HREAL, 1 3X , 4H I MAG , 1 5X4HREA L , 13X,4HIMAG) 

2003 FORMAT (25H AZMUITHAL ORDER MC0MP=,I2} 

2004 FORMAT! 1 5 , 2E 17 . 8 , 2X , 2E 17 . 8 ) 

END 
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Printout D-l (contd) 

(c) Subroutine FIELDS (far-field approximation) 

$ I BF TC FIELD 

SUBROUTINE FIELDS 

THIS SUBROUTINE ASSUMES FAR-FIELD BEHAVIOR FOR THE INCIDENT FIELDS 
TO OBTAIN THE MAGNETIC FIELDS ON A SURFACE S 
A. LUDWIG MARCH 1969 

DIMENSION TITLE{13),PSI(l8l)»E(l8l)tEP{181),H(181),HP{181) 
DIMENSION dUM(2) 

COMPLEX EE»EXPJX 

EQUIVALENCE { , DUM ( 1) ,CO) , (DUM(2) ,SI ) 

COMMON/ GR I D 1 / S I T ( 15) ,COT{ 1 5 I , S I P ( 18 1 ) ,COP ( 1 8 1 ) , T < 15) ,P( 181 ) 
COMMON/GRlD2/SITT( 181) ,COTT( 181 ) ,SIPP { 8 ) ,COPP ( 8 ) , TT( 18 1 } , PP i 8 ) 

DTR=0. 017^53293 
WRITE(6,2001 ) 

01 FORMAT ( 52H I LLUMAT I ON FIELD FROM CIRCULARILY SYMMETRIC FEED//) 

READ IN INPUT FIELDS 
READ(5,1001) TITLE 

READ! 5, 1002) JINM, JO, JIN, IC1 , IC2,MC0MP 
FMC-McOMP 

READ ( 5, 1003) (PSI (J),E(J) ,EP(J) , H ( J ) , HP { J ) , J = 1 , J I N ) 

FOR IC1 LESS THAN OR EQUAL TO 0 CONVERT FROM DB TO VOLTS 
IF(ICl) 10,10,20 
10 DO 15 J = 1 , J I N 

E(J)=10.0**(E< JI/20.0) 

15 H(J)=10.0**(H( J)/20*0) 

20 CONTINUE 

FOR I C 2 GREATER THAN ZERO NEGLECT FEED PHASE PATTERN 
IF { IC2 ) 40,40, 30 
30 DO 35 J= 1 , J I N 
EP( J)=0 
35 HP ( J ) = 0 
40 CONTINUE 

PRINT OUT FIELD ILLUMINATION 
WRITE (6 ,2002)TITLE 

WRITE (6, 2003) ( PS I ( J ) , E ( J ) , EP ( J ) , H ( J ) , HP ( J ) , J =1 , J I N ) 

CONVERT POLAR ANGLES TO RADIANS 
DO 45 J = 1 , J I N 
45 PSI ( J)=DTR*PSI ( J) 

CONVERT TO REAL AND IMAGINARY 
IF { IC2 ) 50 , 50 ,60 
50 DO 55 J= 1 » J I N 
TH=DT R*EP ( J ) 

EE=EXPJX(TH) 

E P ( J ) -E ( J ) *S I 
E( J)-E(J)*CO 
TH=DTR*HP ( J ) 

EE=EXP JX ( TH) 

HP( J)=HLJ)*SI 
55 H(J)=H(J)*CO 
60 CONTINUE 

SET FIELDS OUTSIDE OF INPUT' RANGE EQUAL TO ZERO 
I F ( J I NM- J I N- 1 ) 7 2 , 7 1 » 7 0 

70 PS I ( J IN+2 ) =3. 142 
E ( J I N+2 ) =0 
EP ( J I N+2 ) =0 
H ( J I N+2 ) =0 
HP { J I N+2 ) = 0 

71 PSI ( JIN+1)=PSI ( JINI+.01 
E( J1N+1 )=0 
EP ( J IN + 1 ) =0 
H( JIN + 1 ) =0 
HP ( J I N + 1 ) =0 
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on no o on 


Printout D-l (contd) 


72 CONTINUE 
RETURN 

ENTER HERE FOR GRIDl FIELDS 

ENTRY FIELD1(MMAX,NMAX,HR,HT,HF,R ) 

DIMENSION R ( 15,91) 

COMPLEX HR ( 15,91 ) ,HT( 15,91) ,HF( 15,91 ) 

LINEARIlY INTERPOLATE FIELDS TO INTEGRATION GRID POLAR ANGLES 
J = 1 

DO 200 M=1,mMAX 

109 I F { P S I (J)-T(m) ) 110,110,120 

110 J = J + 1 
GO TO 109 

120 F=(T ( M ) -P $ I (J-l) )/( PSI(J)— PSI ( J - 1 ) ) 

0F=1.0-F 

ETR =F*E( J)+0F*E( J-l) 

ETI =F*EP( J)+0F*EP( J-l) 

.EPR = F*H( J)+OF*H( J-l) 

100 EPI =F*HP ( J ) +0F*HP ( J — 1 ) 

DO 200 N=1,NMAX 
T 1=ETR *SIN{FMC*P(N) ) 

T2=ETI *$ I N ( FMC *P ( N ) ) 

T3=EP R *COS(FMC*P(N) ) 

T4= EP I *CO$(FMC*P(N)) 

ASSUME FAR FIELD BEHAVIOR HR=0 , HP= ET , HT=-E P 
HR(M,N)=(0. 0,0.0) 

HT(M,N)=CMPLX(-T3,-T4) 

HF(M,N)=CMPLX( T1,T2) 

200 CONTINUE 
RETURN 

ENTER HERE FOR GRID2 FIELDS 
ENTRY FIELD2(JQ»K0»ETT0,EPPQ) 

COMPLEX ETTO , EPPO 

LINEARLY INTERPOLATE FIELDS TO OUTPUT GRID POLAR ANGLES 
J=1 
J J = JO 

509 I F ( PS I ( J ) -TT { J J ) )510,510,520 

510 I F ( J • EQ • 18 1 ) GO TO 520 
J=J + 1 
GO TO 509 

520 F=(TT( JJ)“PSI ( J-l ) ) /(PSH J)-PSI ( J-l) ) 

0F=1 . 0- F 

ETR = F*E( J)+OF*E( J-l) 

ETI =F*EP( J)+QF*EP< J-l ) 

EPR' = F*H( J)+GF*H( J-l) 

500 EPI = F*HP { J ) +OF#HP { J-l ) 

C 

K=KO 

Tl= ETR *SIN(FMC*PP(K) ) 

T2 = ETI * S I N ( F MC * P P ( K ) ) 

T3=EPR *COS( FMC*PP{K) ) 

T A=EP I *COS(FMC*PP(K) ) 

ETT0=CMPLX(T1,T2) 

EPPO=CMPLX (T3,T4) 

RETURN 

1001 FORMAT { 1 3A6 ) 

1002 FORMAT (1015) 

1003 FORMAT ( 5 F 1 0 • 6 ) 

2002 FORMAT ( 13A6/47H POLAR E-PLANE H-PLANE/ 

150H ANGLE VOLTS DEG VOLTS DEG) 

2003 FORMAT (F10.2,F12.6,F8.2,F13.6,F8.2) 

END 
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Printout D-l (contd) 
(d) Subroutine SURF 


$ IBFTC SURFDK 

subroutine SURf(MMA X »NMAX,F,FT,FP) 

C 

C THIS SUBROUTINE PROVIDES TrE MAIN PROGRAM WITH RHO AND ITS 

G PARTIAL DERIVATIVES WITH RESPECT TO THETA AND PHI 

C RHO IS DENOTED BY p IN THIS PROGRAM FOR SOME OBSCURE REASON 

C 

DIMENSION TIN! 15) ,A0(15) ,A1( 15) ,A2( 15) , DAO ( 1 5 ) , DAI { 1 5 ) ,DA2( 15) 
DIMENSION TITLE! 12) 

COMMON/GRId 1/SIT( 15) ,COT! 15) ,Sj P ( 181 ) f COP( 181 ) ,T< 15 ) , P <18 1 ) 
DIMENSION F(1&,91 ),fT( 15,91) ,FP! 15,91) 

C 

C READ IN FOURIER C OE I FF IC I ENTS 

C 

4 READ ! 5 » 100DTITLE 
WRI T fc ( 6 , 2001ITITLE 

READ (5, 1002) (TIN(M) , AO(M) ,A1(M) ,A2(M) ,DAO(M) ,DA1(M) ,UA2(M) , 

*M= It MM AX) 

WRITE (6 , 2003) ( T I N ( M ) , AO ( M ) , A1 ( M ) , A2 ( M ) , DAO ! M ) , DAI ( M ) , DA2 i M ) , 
*M=1,MMAX) 

C 

C CHECK THAT INPUT DATA AND INTEGRATION GRID AGREE 

C 

DO 10 M= 1 f MM AX 
TCHE=T (M)/0*017453293 

IF! ABS(TCHE-TIN(M) ) 00000 1 ) 10 , 99 , 99 
99 WRITE (6 ,2002) TIN (M) , TCHE 
10 CONTINUE 
C 

C FILL OUT SURFACE TABLE 

C 

DO 50 N=1»NMAX 
DO 50 M=1,MMAX 
C0P2P=1*0-2.0*SIP(N)*SIP(N) 

F(M,N)=A0(M)+A1(M)*C0P(N)+A2(M)*C0P2P 
FT(M,N)=DA0<M)+DA1 (M) *C0p! N)+DA2 (M)*C0P2P 
SIP2P=2.0 *SIP(n)*CgP(N) 

FP(M,N)= -A1(m>*SIP(N)-2.0*A2(M)*SIP2P 
50 CONTINUE 
RETURN 

1001 FORMAT ! 1 2 A6 ) 

1002 FORMAT! 7F10. 4) 

2001 FORMAT ! 44H0 SURFACE DAT. A INPUT IN THE FORM OF FOURIER/ 

*37H COEFFICIENTS AND THEIR DER I VAT I VES//12A6// 

*44H THETA AO ! THETA ) A1 ! THETA) A2 (THETA ) » 

*37H D AO/DTHETA DA1/DTHETA DA2/DTHET A ) 

2002 FORMAT (55H SURFACE INPUT DATA INCONSISTANT WITH INTEGRATION GRID/ 
*14H INPUT ANGLE, F7. 2, 18H DEG GRID ANGL E , F7 • 2 , 4H DEG) 

2003 FORMAT (F8*2,6F12. 4) 

END 


(e) Subroutine SETUP 


$ IBFTC SETO 

subroutine setup (ngi, i, jmax,kmax,mmax, nmax> 


c 

c 

c 

c 


THIS PROGRAM ESTABLISHES ThE OUTPUT GRID 
GRIDS, AND PRECOMPUTES TRIG FUNCTIONS ON 


AND UP TO 5 INTEGRATION 
ALL OF THE GRID POINTS 


DIMENSION MM! 5 ) , NN ( 5 ) , T! 15, 5),P(181, 5) 

COMMON /GR ID 1/S I T ( 15), COT! 15) , SIP! 181 ), COP! 181) , TGI ( 15 ) , PG 1 ( 181) 

COMMON/GRI D2/SI TT ! 18 1 ) ,COTT( 181 ) ,SIPP( 8) ,COPP( 8 ) , TT ( 181 ), PP !8) 
DIMENSION dUM(2) 

COMPLEX E , EXP JX 

EQUIVALENCE ( E , DUM ! 1 ) ,C0), (DUM(2) ,SI) 

DATA LI / 06 06 060606374/, L2/ 06060 606 04 77 4/, L3/ 06 0606063 63 74/ 

DATA L4/ 0606060474774/ , L 5/0341 3qq000000/ 
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Printout D-l (contd) 


C GRID VALUES MAY BE READ IN BY INPUTING ONLY -MAX VALUES, OR IF 

C EQUALLY SPACED VALUES ARE DESIRED, PROGRAM WILL COMPUTE VALUES 

C GIVEN A STARTING POINT AND AN INCREMENT 

C 
C 

C READ IN OUTPUT GRID 

READ(5,1001)JMAX»TT1» DTT 
IF(DTT) 10,20,10 
10 DO 15 J= 1 , JMAX 
F J= J— 1 

15 TT(J)=TT1+FJ*DTT 
GO TO 30 

20 READ(5,1002)(TT(J),J=1,JMAX) 

30 RE AD (5, 1001)KMAX,PP1»DPP 
I F ( DPP ) 40 ,50,40 
40 DO 45 K= 1 , KMAX 
FK=K-1 

45 PP(K)=PP1+FK*DPP 
GO TO 60 

50 READ(5, 1002) (PP(K) ,K=1, KMAX) 

60 CONTINUE 
C 

C read in integration GRIDS 

READ ( 5 , 1001 ) NG1 
DO 100 1=1, NG1 
RE AD ( 5 , 1001 ) MM ( I ),T1,DT 
MMM=MM ( I ) 

IF(DT) 110,120, 110 
110 DO 115 M= 1 , MMM 
FM=M-1 

1 1 5 T ( M , I )=T1+FM*DT 
GO TO 130 

120 RE AD (5, 1002) ( T ( M , I ),M=1,MMM) 

130 READ(5,1001)NN( I ) ,Pl,Dp 
NNN=NN ( I ) 

IF(DP) 140,150,140 
140 DO 145 N=1 , NNN 
FN=N-1 

145 P ( N , I ) = P1 + FN*DP 
GO TO 100 

150 RE AD (5, 1002) ( P ( N , I ),N=1»NNN) 

100 CONTINUE 
C 

C PRINT OUT GRID DATA 

WR I TE (6 , 200 1 ) 

WRITE (6, 2006) ( L3 , J , L5 , TT ( J ) , J= 1 , JMAX ) 

WR I TE ( 6 , 2002 ) 

WRITE (6, 2006) ( L4, K , L5 , PP ( K ) , K= 1 , KMAX ) 

WRITE(6,2004)NG1 
DO 155 1=1, NG1 
WRITE(6, 2005)1 
MMM=MM ( I ) 

NNN=NN ( I ) 

WRITE (6, 2006) ( L1,M,L5,T(M,I) ,M=1,MMM) 

WR I TE ( 6 , 2002 ) 

155 WR ITE (6,2006 ) (L2,N,L5,P(N, I ) ,N=1,NNN) 

C 

C establish GRI D2 TABLE 
DTR=0. 017453293 
DO 200 J=1,JMAX 
TT( J)=DTR*TT( J) 

E = EXP JX ( TT ( J ) ) 

S I TT ( J ) = S I 
200 COTT ( J ) =C0 

DO 210 K= 1 , KMAX 
PP(K)=DTR*PP(K) 

E=EXPJX ( PP ( K ) ) 

S I PP ( K ) =S I 
210 COPP ( K )=C0 
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Printout D~1 (contd) 


c establish gridi tables 
i=o 
c 
c 

C FOR NG1 GREATER THAN I* RE-ENTEr SUBROUTINE HERE FOR NEW GRlDl 

C 

ENTRY RESET( I , MMAX, NMAX) 

1 = 1 + 1 

MMAX=MM ( I ) 

DO 400 M= 1 » MMAX 
TG1(M)=DTR*T(M,I ) 

E=EXP JX ( TG 1 { M ) ) 

SIT ( M ) =S I 
400 COT ( M ) =C0 
NMAX=NN ( I ) 

DO 410 N=1,NMAX 
PG1 (N)=DTR*P(N, I ) 

E = EXP JX { PG 1 ( N ) ) 

SIP (N)=SI 
410 COP ( N ) =C0 
RETURN 

1001 FORMAT ( I5»2F10.0) 

1002 F0RMAT(8F10.0) 

2001 FORMAT (50H THE FOLLOWING OUTPUT GRID HAS BEEN ESTABLISHED/ 

125H ALL ANGLES IN DEGREES//) 

2002 FORMAT ( // ) 

2004 FORMAT { 17H0 THE FOL L OW I NG , I 3 , 40H INTEGRATION GRIDS HAVE BEEN EST 
1ABLISHED/) 

2005 FORMAT ( 1 8H0 SEGMENT NUMBER ,13/) 

2006 FORMAT (5(A6,I3,A2>F9*4) ) 

END 


(f) Subroutine PATHL 


$ I B FTC PATHLD 

SUBROUTINE PAThL {RHO, J,K, MMAX, NMAX, GAM) 

C THIS SUBROUTINE EVALUATES ThE PATH LENGTH FUNCTION GaMMA 

DIMENSION RHO ( 15, 91),GAM( 15, 91) 

COMMON/ GR I Dl/SI T ( 15) ,COT( 15 ) , SI P ( 181 ) ,CUP { 181 ) 

COMMON /GRID 2/ SI TT ( 181 ) ,COTT ( 181 ) , SI PP ( 8) ,COPP( 8 ) » TTU81) «PP <8 ) 

DO 10 M-l , MMAX 

Tl=SIj (M)*SITT( J )*CoPP(K) 

T2=SIT(M)*SITT( J)*SIPP(K) 

T3=C0t<M)*C0TT( JJ-1.0 
DO 10 N = 1 , NMAX 

10 GAM(M,N)=RH0(M f N>*(Tl*C0P{N)+T2*SIP(N)+T3) 

RETURN 

END 


(g) Subroutine FINT 


$ I BFTC FINtD 

SUBROUTINE F I NT ( X , Y , f , R , MMAX , NMAX , STOT ) 

C 

C THIS SUBROUTINE NUMERICALLY INTEGRATES 

C THE RADIATION INTEGRAL 

C A* LUDWIG JUNE 1968 

C 

DIMENSION X< 1 5 ) t Y ( 18 1 ) ,R( 15, 9 1 ) 

DIMENSION DUM( 2) 

COMPLEX F ( 15, 91,3) , STOT (3) ,E,T1,T2,T3,A,B,C,EXPJX, F12,F23, 

1F14 , SUM, TOT 

EQUIVALENCE ( E , DUM ( 1 ) , CO ) , ( DUM ( 2 ) , S I ) 

DO 10 L= 1 , 3 
10 ST0T(L)={0. 0,0.0) 

DO 200 N= 1 » NMA X 
DY=0.5*(Y(N+1)-Y(N) ) 
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Printout D-l (contd) 


DO 200 M- 1 ,MMAX 
DS=DY*{X(M+1)-X{M) ) 

R 1 = R { M + l , N+ 1 ) ~R ( M , N ) 

R2=R(M,N+1)-R(M+1,N) 

R3=R(M,N)+R(M+1,N) 

BE= ( R 1 + R2 ) *0 * 5 
CE = ( R 1-*R 2 ) *0 • 5 
AL= ( R3-C E ) *0 • 5 
I F ( ABS( BE >-0.01) 100,100, 110 
100 F3I=BE*0. 33333333 
F1I=BE*0.5 
F1R=1.0-F1I*F3I 
F3R=0.5-BE*BE/8.0 
GO TO 140 
110 E = EXP J X ( BE ) 

F 1R=S I / BE 
FllMl.O-COJ/BE 
F3R = F 1R-F 1 1 /BE 
F3I=<F1R-CQ)/BE 

140 IF( ABS(CE)-O.Ol) 150,150, 160 
150 F4I=CE*0 .33333333 
F 2 1 =CE*0 • 5 
F2R=i.0-F2I*F4l 
F4R=0.5-CE*CE/8.0 
GO TO 170 
160 E = EXPJX ( C E) 

F 2 R=SI/CE 
F2I=(1.0-C0)/Ce 
F4R=F2R-F2 I /C E 
F4I = ( F 2R-C0 ) /C E 
170 E=EXPJX(AL) 

F12=CMPlX(F1R*F2R-F1I*F2I , FI I*F2R+F1R*F2 I ) 
F23=CMPLX(F2R*F3R-F2I*F3I , F2 I*F3r+F2R*F3 I ) 
F14=CMPlX{F1R*F4R-f1I*F4I,F1I*F4R+F1R*F4I ) 
DO 200 L = 1 , 3 
T1=F{M+1,N+1 ,L)-F(M,N,L) 

T 2 =F (M,N+1,L)-F(M+1,N,L) 
T3=F(M,N,L)+F(M+1,N,L) 

B=T1+T2 
C=T 1-T 2 
A=T3-0.5*C 

$UM=A*F12+B*F23+C*Fl4 
TOT=E*SUM*DS 
200 STOT ( L )=STQT ( L ) +TOT 
RETURN 
END 


(h) Subroutine PRTIM 


$ I BFTC PRTM13 PRINT TIMERS 

SUBROUTINE PRTImI 

PRTM 

20 

C VERSION 2, FEBRUARY 1, 1967 

PRTM 

30 

C VERSION 3, APRIL 1968. 

TOT AL=0. 

PRTM 

40 

PART =0 

PRTM 

50 

CALL GETTIM (NTIME) 

PRTM 

60 

GO TO 10 

PRTM 

70 

C 

PRTM 

80 

c 

PRTM 

90 

c 

PRTM 

100 

ENTRY PRTIM2 

PRTM 

110 

CALL GETTIm(ITIME) 

PRTM 

120 

T0TAL=FL0ATUTIMe-NTIME)/3600. 

PRTM 

130 

PART = FlOAT( I TIME- 1 PART)/3600. 

PRTM 

140 

10 TOTSEC = T0TAL*60 • 

PARSEC = PART*60. 

WR I TE{ 6, 1000) PART, PARS EC , TOTAL, TOTSEC 

CALL GETTIm( IPART) 

PRTM 

160 

RETURN 

PRTM 

170 

1000 FORMAT ( 15H0PARTIAL TIME =Fl0.4,5H M I NS , 3X , 1H= , F 1 2 .4 , 5H SEqS, 
* 10X , 12HT0T AL TIME =F10.4,5h MI NS , 3X , 1H= , F12 .4, 5H SECS) 

END 

PRTM 

200 


50 
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Printout D-l (contd) 


$ I BMA P 

FSTIME 


FETCH AND STORE LOC 5 

FST I 

10 

* 

FSTIME 


FETCH AND STORE LOC 5 

FST I 

20 





FST I 

30 

* 



CALL SETTIM{ IT IMF ) 

FSTI 

40 

* 



CALL GETTIM! ITIME) 

FST I 

50 

T 




FSTI 

60 


ENTRY 

SETTIM 


FSTI 

70 


ENTRY 

GETTIM 


FSTI 

80 

SETTIM 

CLA* 

3,4 


FSTI 

90 


STO 

5 


FSTI 

100 


TRA 

1,4 


FSTI 

110 

GETTIM 

CLA 

5 


FSTI 

120 


S TO* 

3,4 


FSTI 

130 


TRA 

1,4 


FSTI 

140 


END 



FSTI 

150 
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Printout D-2. Computer printout of spherical-wave expansion program 
(a) Main program (orthogonality technique) 


$ I BFTC MAINDK 

C SPHERICAL WAVE EXPANSION PROGRAM - ORTHOGONALITY TECHNIQUE 

c this program determines the spherical wave coefficients for the 

C M-TH FOURIER COMPONENT OF A RADIATION PATTERN. (SEE STRATTON P 416) 

C COEFFICIENTS ARE OBTAINED FROM ORTHOGONALITY INTEGRALS 

C THE EXPANSION MATCHES THE FAR-F J ELD INPUT PATTERN WITH THE FAR- 

C FIELD FORM OF SPHERICAL WAVES. ART LUDWIG 12/23/68 

C 


DIMENSION NAME ( 13) , NAME J { 13) ,NDUM{ 13) 
DIMENSION PA(60i ,PB(60) 

DIMENSION T( 181 ) t E < 181 ) i EP ( 181 ) ,H(181 
DIMENSION F (60* 18 1 ) » G ( 60 i 1 8 1 ) , PM I 6 1 ) 
DIMENSION A( 181 9 2 ) 9 B ( 18 1 1 2 ) 

DIMENSION ACOEI60,2) ,BC0E(60,2) 
DIMENSION AOUT I 1 , 2 ) , BOUT 11,2) 

DATA NDUM/ 1 3^0606060606060/ 

DTR=0. 017453293 


» HP ( 18 1 ) 


C 

C 

C 


C 

C 


C 


C 


C 


C 

C 

C 


READ IN M-TH COMPONENT OF INPUT PATTERN 

1 READI5,1001)NAMEJ 
CALL PRTIM1 

READ(5,1002)MCOMP,NMAX 

READ (5, 1001) NAME 

READ (5, 1002) JMAX,J 0 ,JIN, IC1, IC2 

READ (5, 1003) (T(M) ,E(M),EP(M>,HIM),HPIM),M=1,JIN} 

FOR IC1 LESS THAN OR EQUAL TO 0 CONVERT FROM DB TO VOLTS 
IF ( I Cl ) 10 » 10 9 20 
10 DO 15 J = 1 , J I N 

E( J)=10.0**(E< J)/ 20 . 0 ) 

15 HI J)=10.0**(H( JJ/20.0) 

20 CONTINUE 

FOR IC2 GREATER THAN 0 NEGLECT PHASE 
IF (I C2) 40, 40, 30 
30 DO 35 J=1,JIN 
EP( j)= 0 . 

35 HP ( J ) =0. 

40 CONTINUE 

PRINT OUT INPUT PATTERN 
WRITE(6,2001)NAMEJ 
WR I TE ( 6 , 2006 ) MC QMP 
WRITE(6,2002)NAME 

WRITE (6, 2003) (TIM) ,E(M),EPIM),H(M) , HP I M) ,M = 1 , J IN ) 

CONVERT TO RADIANS AND REAL AND IMAGINARY 
DO 45 J=l, JIN 
45 T( J)=DtR*T( J) 

IF I I C2 ) 50 , 50 , 60 
50 DO 55 J= 1 , J I N 
TH=DtR*EP I J ) 

EP I J)=£( JFSINITH) 

El JFEI JFCOSITH) 

TH=DTR*HP{ J) 

HP(J)=HIJ)*SIN(TH) 

55 HI J)=H( JFCOSITH) 

60 CONTINUE 

Establish ai theta) dtheta and bithetaidtheta vectors 

All, 1)=E( 1)*T( 1) 

A( It 2)*EP( 1)*T( 1) 

B I 1, 1)=HI 1)*T( 1) 

BU,2) = HP( 1)*T( 1) 

DO 65 J=2,JIN 
DTH=T(J)-T( J-l) 

A( J, 1)=E( J)*DTH 
A( J,2) = EP( JFDTH 
B( J, 1)=H( J)*DTH 
65 B( J,2)=HPI JFDTH 
CALL PRTIM2 
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Printout D-2 (contd) 


OBTAIN F AND G MATRICES 

FMC^MCOMP 
DO 80 J= 1 » J I N 
Z=C0S I T I J ) ) 

DO 83 N=l, MCOMP 
85 PM IN ) =0 ♦ 

NC=NMAX+1 

CALL LEGEND INC , MC OMP , Z , PM ) 

DO 80 I = 1 » NM AX 
FII , J ) = PM ( I )*FMC 
T 1= I -MCOMP+1 

T2=I+1 

G ( I , .J)«TI*Pm( I + l)-T2*Z*pM( I ) 

80 CONTINUE 

PERFORM NUMERICAL INTEGRATION BY MATRIX MULTIPLICATION 

CALL MULTI JIN, NMAX,2»F» A, ACOE, 0,60, 181 , 181,2,60,2) 

CALL MULTI JIN, NMAX, 2, G,B, ACOE, 1,60, 181,181,2,60,2) 

CALL MULTI JIN,NMAX, 2, F,B,BCOE, 0,60, 181, 181,2,60,2) 

CALL MULTI JIN, NMAX,2,G, A, BCOE, 1,60, 181,181,2,60,2) 

NORMALIZE C OE I FF IC I ENTS AND COMPUTE POWERS 

P 10 2Z= (3. 1415 92 7/2. 0)*0 .002 655 
PT0T = 0 * 

DO 95 N- 1 , NM AX 
P A IN ) =0 • 

P B I N ) ~ 0 ♦ 

F ACT = 1 * 

I F I MCOMP ) 92 , 9 2 , 90 

90 DO 91 M=l, MCOMP 
FF=MN-M+1)*(N+M) 

91 FACT=FACT*FF 

92 FF = 2'MM+1 

fact=fact/ff 

FF=2*N*(N+1) 

FACT=FACT*FF 
K = 0 

93 K=K+1 

BT^BCUE IN,K) 

AT=ACUE I N, K i 

BCOE(N,K)=BCOEIN,K)/FACT 
ACOE I N ,K ) = ACOE I N t K ) / FACT 
PA(N)=PA(N )+AT^ACOElN, K)*pI02Z 
PB(N}=PBIN)+BT*BC0E{N,K)*PI02Z 
IF! K- 1)93, 93, 94 

94 CONTINUE 
PTOT=PTOT+PA(N)+PB(N) 

95 CONTINUE 

OUTPUT C OE I FF IC I ENTS 

WRI TEI6, 2001 )NAMEJ 
PUNCH 1001 , NAME J 
PUNCH 1001, NAME 
PUNCH 2010 ,NMAX, MCOMP 
WRI TE ( 6, 2004) MCOMP 
P$UM=0 

DO 100 J= 1 , NMAX 
PAIJ)=PAIJ ) /PTOT 
PBI J)=PB( J)/PTOT 
PSUM = PSUM+P A I J ) +PB I J ) 

PUNCH 2005, J , ACOE I J , 1 ) , AC OE I J , 2 ) , BC OE ( J , 1 ) , BC OE I J , 2 ) 

100 WRITE 1 6, 2005) J,ACOE( J, 1) ,ACOEl J,2) ,BCOE( J,l) ,BCOEI J,2),PA( J ) , PB I J) 
* , P SUM 

WRITE (6,2007)PTOT 
CALL PRTIM2 



O o o 


Printout D-2 (contd) 


OUTPUT SUMMATION OF SPHERICAL MODES 

REA0(5, 1002) JMAX, JIN 
J0 = 0 
IC1=1 
I C2 = -l 

PUNCH 1001,NAMEJ 

PUNCH 1002, JMAX, JO, JIN, IC1, IC2 

WRITE (6, 2008 )NAMEJ 

WRI TE( 6, 2002 )N0UM 

DT= JMAX-1 

DT= 180 *0/DT 

J = 1 

TH = 0 

I F { MCOM P-1 ) 170, 180,170 
170 DO 175 N=1,NMAX 
F ( 1 , N ) = 0 
175 G(1,N)=0 
GO TO 215 

180 DO 190 N= 1 , NMAX 
FN=N*(N+1) 

F ( 1 >N) -FN/2.0 
190 G(l,N)=FN/2.0 
GO TO 215 
200 FJ=J-1 
TH=F J*DT 
Z=COS { TH*DTR ) 

S=S IN( TH^DTR ) 

DO 205 N*1,MC0MP 
205 PM { N ) =0 

CALL LEGEND { NC » MC OMP , Z , PM ) 

DO 210 N- 1 » NMAX 
F(1,N)=FMC*PM(N)/S 
T1=N“MC0MP+1 
T2=N+1 

G{ i,N) = t1* pm ( n+ D-T2-Z-PM( N) 

210 G(1,N)=G( 1,N)/S 

215 CALL MULT{NMAX,1,2,F,AC0E,A0UT,0, 60,181, 60,2,1,2) 

CALL MUlT(NmAX,1,2,G,BC0E,A0UT,1, 60,181, 60,2,1,2) 

CALL MULT ( NMAX , 1 , 2 , G , AC OE , BOUT , 0 , 60,181, 60,2,1,2) 

CALL MUlT ( Nm AX , l , 2 , F , BC OE , BOUT , 1 , 60,181, 60,2,1,2) 

CALL VECTOR (AOUTt 1, 1 ) ,AOUT( 1,2) ,EAMP,EPHI ) 

CALL VECTOR ( BOUT (1,1 ) ,BOUT( 1 ,2) ,HAMP,HPHI ) 

PUNCH 2009, TH,EAMP»EPHI ,HAMP,HPHI 
WRITE(6,2003) TH,EAMP,EPHI ,HAMP,HPHI 
J - J+ 1 

I F ( J- J I N ) 200 ,280,300 
280 IF { J- JMAX ) 200 , 285 , 300 
285 I S I GN= 1 
TH=180* 

IF (MCOMP-1) 370,380,370 
370 DO 375 N=1,NMAX 
F ( 1 ,N ) =0 
375 G ( 1 , N ) =0 
GO TU 215 

380 DO 290 N= 1 , NMAX 

isign=-isign 

FN=N S MN + 1 ) * I S IGN 
F(l,N)=-FN/2.0 
290 G ( l,N)=FN/2.0 
GO TO 215 
300 CONTINUE 

CALL PRTIM2 
GO TO 1 

1001 FORMAT (13A6) 

1002 FORMAT (1015) 

1003 FORMAT (5F10.0) 

2001 FORMAT { 1H1,13A6) 

2002 FORMAT (1H0, 13A6/46H POLAR E-PLANE H-PLANE/ 
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1 50H ANGLE VOLTS DEG VOLTS DEG) 

2003 FORMAT (F 10. 2, FI 2.6* F8.2,F13.6,F8.2) 

2004 FORMAT ( 5 1H0 SPHERICAL WAVE COEFFICIENTS FOR AZMUITHAL ORDER, 12// 
120X,4HA(N) ,32X,4HB(N) , 27X , 28HFRACT ION OF TOTAL MODE POWER/ 

25H N, 7X, 4HREAL, 13x,4HIMAG, I5X,4HREAL, 13X,4HIMAG, 

3 I 3X , 7HA MODES, 6X,7HB MODES , 5X , 16HCUMULAT I VE TOTAL) 

2005 FORMAT (I5,2E17*8,2X,2E17.8,2X,2E14.5,F14.8) 

2006 FORMAT { 49H0 INPUT PATTERN FOR AZMUITHAL COMPONENT OF ORDER, 12) 

2007 FORMAt(1HO,I9H ToTAL MODE POWER , E 15 . 8 , 6H WATTS) 

2008 FORMAT ( 1H1, 13A6/41H FAR FIELD SUMMATION OF SPHERICAL MODES) 

2009 FORMAT (F 10. 2 , F 10. 6 , F10. 2 , F10 . 6 , F10. 2 > 

2010 FORMAT (215) 

END 


(b) Main program (linear equation technique) 

SIBFTC MAINDK 

c spherical wave e xpansion program - linear equation technique 

THIS PROGRAM Determines ThE SPHERICAL WAVE C OE I FF 1C I ENTS FOR THE 
M-TH FOURIER COMPONENT OF A RADIATION PATTERN. (SEE STRATTON P 416) 
THE COEFFICIENTS ARE OBTAINED BY INVERTING A SYSTEM OF LINEAR 
EQUATIONS- WARNING- THE INPUT DATA POINTS SHOULD BE NEARLY EQUALLY 
SPACED BETWEEN 0 AND 180 DEGREES. OTHERWISE THE PROBLEM BECOMES 
VERY ILL CONDITIONED 

THE EXPANSION MATCHES THE FAR-F I ELD INPUT PATTERN WITH THE FAR- 
FIELD FORM OF SPHERICAL WAVES. ArT LUDWIG 12/23/68 

DIMENSION NAME ( 13 ) 

DIMENSION T( 181) ,E{ 181) ,EP( 181 ) ,H( 181 ) ,HP( 181) 

DIMENSION A ( 100 , 100 ) ,XR( 100) ,X I ( loO) »BR( 100) ,BI ( 100) , PM (50) 
DIMENSION C ( 100 , 100 ) , S ( 500 ) 

DTR=0. 017453293 

READ IN M-TH COMPONENT OF INPUT PATTERN 

1 READ ( 5 , 100DNAME 
CALL prtimi 

READ (5, 1002) JMAX, JO, JIN, IC1 , IC2, MCOMP 
RE AD (5, 1003) (T(M) ,E(M) , £P ( M) ,H(M) ,HP(M) ,M=1,JIN) 

FOR IC1 LESS THAN OR EQUAL TO 0 CONVERT FROM DB TO VOLTS 
IF(ICl) 10, 10,20 
10 DO 15 J= 1 , J I N 

E( J)=10.0**(E{ JJ/20.0) 

15 H( J) = 10.0**(HU>/20.0) 

20 CONTINUE 
C FOR IC2 GREATER THAN 0 NEGLECT PHASE 

IF ( IC2)40»40,30 
30 DO 35 J= 1 , J I N 
EP ( J ) =0 • 

35 HP ( J )=0. 

40 CONTINUE 

C PRINT OUT INPUT PATTERN 

WRITE(6,2002)NAME 

WRITE (6, 2003) (T(M),E(M),EP(M),H(M),HP(M),M=1,JIN) 

C CONVERT TO RADIANS AND REAL AND IMAGINARY 

00 45 J= 1 , J I N 
45 T( J)=DTR*T( J) 

IF ( IC2) 50, 50 ,60 
50 DO 55 J= 1 , J I N 
TH=DTR*EP ( J ) 

EP( J)=E( J)*SIN(TH) 

E( J)=E( J)«COS(TH) 

TH=DTR*HP( J) 

HP(J)=H( J)*SIN(TH) 

55 H( J)=H(J)*COS(TH) 

60 CONTINUE 
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c 

C THE COE IFF I C I ENT VECTOR X IS OBTAINED FROM THE MATRIX EQUATION 

C AX=B, WHERE THE B VECTOR IS DEFINED IN THE FOLLOWING SECTION 

C BR AND B I ARE REAL AND IMAGINARY PARTS OF B 

C 

DO TO J=1,JIN 
BR(2*J~1)=E( J) 

BR!2*J)=H( J) 

81 ( 2*J-1 ) = EP! J ) 

70 B I ( 2* J ) =HP { J ) 

C 

C THE ELEMENTS OF THE MATRIX A ARE OBTAINED FROM ASSOCIATED 

C LEGENDRE POLYNOMIALS AND THEIR DERIVATIVES AS FOLLOWS 

C 

CALL PRTIM2 
FMC=MCOMP 
DO 80 K-ltJIN 
Z=COS { T ( K ) I 
S=$ IN ( T( K ) ) 

DO 85 N=1,MC0MP 
85 PM (N ) -0 • 

NMAX=JIN+1 

CALL LEGEND ( NMAX , MCOMP » Z » PM ) 

DO 80 J= 1 , J I N 
F=FMC*PM( J)/S 
Tl- J-McOMP+1 
T2=J+1 

G=T1*PM! J+1)-T2*Z*PM( J) 

G=G/S 

A(2*K-1,2*J-1)=F 
A(2*Kr2*J)=F 
A ( 2*K- 1 » 2* J ) =G 
A(2*K»2*J-1)=G 
80 CONTINUE 
C 

C THE SYSTEM IS INVERTED TO YIELD X 

C 

EPS=. 00000001 
NS0LVE=2*JIN 
CALL PRTIM2 

CALL SOLVE! 100 , A , NSOLV E , BR , XR , 1 , EPS , 2 , IT1,C ,S ) 

CALL SOLVE! 1 00 , A , NSoLVE , B I ,XI ,2, EPS, 2, IT2,C,S) 

CALL PRTIM2 
IF! I T 1 ) 99 » 99 , 90 
90 IF(IT2)99,99,95 
99 WR I TE ( 6 , 2010 ) IT1,IT2 
95 CONTINUE 
C 

C OUTPUT RESULTS 

C 

PUNCH 1001, NAME 
WRI TE! 6, 2004) 

DO 100 J=1 , J IN 

PUNCH 2005, J,XR(2*J-1) , X I ( 2* J-l ) , XR ( 2* J ) , XI ( 2* J ) 

100 WRI TE<6,2005) J,XR(2*J-1) ,XI! 2* J-l ) ,XR!2*J) ,XI (2*J) 

GO TO 1 

1001 FORMAT ! 1 3 A6 ) 

1002 FORMAT! 1015) 

1003 FORMAT (5F10.0) 

2002 FORMAT! 1H1, 13A6/46H POLAR E-PLANE H-PLANE/ 

1 50H ANGLE VOLTS DEG VOLTS DEG) 

2003 FORMAT! F10 . 2 , FI 2. 6 » F8 . 2 , F 1 3 . 6 , F 8 . 2 ) 

2004 FORMAT (1H1,29H SPHERICAL WAVE C OE I FF I ENTS// 

120X ,4HA ( N ) , 26X , 4HB ( N ) / 

25H N,7X,4HRFAL, 13x»4HIMAG, 15X,4HREAL, 13X,4HIMAG/) 

2005 FORMAT! I 5 , 2E 1 7 . 8 , 2X , 2E 17 ♦ 8 ) 

2010 FORMAT <47H CAUTION-FULL ACCURACY NOT ACHIEVED IN MATRIX/ 

120H INVERSION. IT1=,I2,7H IT2=,I2) 

END 
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(c) Subroutine LEGEND 


$ I BFTC LEGUK FUL 1ST, DECK 

SUBROUTINE LEGENO(NMAX,MfZ»VAL) 

C 

C THIS SUBROUTINE RETURNS VALUES OF THE ASSOCIATED LEGENDRE FUNCTION 

C WITH INTEGER INDICES FROM N=M TO N=NMAX. 

C DO NOT USE NMAX LESS THAN M+l 

C VALUES ARE OBTAINED USING UPWARD RECURSION, AND HAVE CHECKED WITH 

C TABULATED VALUES TO 5 PLACES THRU N-56 FOR M= 1 , AND N=10 FOR H=5 

C HIGHER INDICES WERE NOT CHECKED, 

' C ART LUDWIG 6/6/64 

DIMENSION VALt IOO) 

DOUBLE PRECISION TERM1 , TERM2 , TERM3 , ZD 
ZD = Z 
FM = M 

TERMING. 

I F { M ) 8 , 8 , 9 

8 TERM2= 1,0 
GO TO 11 

9 KMAX=2*M~1 

TERM2- ( 1.0-ZD*ZD)**( FM/2.0) 

DO 10 K= 1 , KM AX » 2 
FK-K 

10 TERM2=TERM2*< 2.0*FM-FK) 

VAL(M)=TERM2 

11 NN=NHAX-1 

DO 20 N=Mf NN 
FN=N 

TERM3-( (2.0*FN+1.0)*ZD*TERM2-( FN+FM ) *TERM1 ) /( FN-FM+1.0) 

VAL ( N+ 1 )- TERM3 
TERM1=TERM2 
20 TERM2=TERM3 
RETURN 
ENO 


(d) Subroutine SPHANK 


SIBFTC SPHAND 

SUBROUTINE SPHANK ( N» R , X *Y ) 

C 

C THE VALUE OF THE SPHERICAL HANKEL FUNCTION IS RETURNED, MULTIPLIED 

C BY THE FACTOR ( - J ) ** ( N+l ) * RHO * EXP(J*RHOJ 

C 

DOUBLE PRECISION TERM, AR , AI 
ar-o 

A I -0 

P 1 = 3* 1415927 
K=0 

TERM= 1 
GO TO 100 
20 K=K+1 
T1=N+K 
T2=N-K+1 
T3=2*K 

TERMITE RM* T1*T2/T3 
TERM-T ERM/R 
GO TO (200, 100) , I GO 
100 AR= AR+TERM 
I G0 = 1 

IF(K-N)20, 1000, 1000 
200 AI =AI -T ERM 
I GO-2 

TERM=-*TERM 
IF (K-N) 20, 1000, 1000 
1000 X = AR 
Y = A I 
RETURN 
END 
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(e) Subroutine VECTOR 

$ I BFTC VDECK 

SUBRQU T INE VEcTOR(X,Y f AMp»PhI > 

this subroutine converts complex data to polar form 

c=o. 

I F { X ) 100,200,300 
lOO IF (Y) 110, 120,120 
HO C = 360 • 

l 20 PHI=A T AN( Y/X)*57.295 7 7951+180.-C 
GO TO 400 

200 IF(Y)210, 220,230 
2 1 0 PH I = — 90 • 

GO TO 400 
220 PH I = 0 • 

GO TO 400 
230 PHI = 90 • 

GO TO 400 

300 PHI=ATAN(Y/X)*57. 29577951 
400 AMP=SQRT( X*X+Y*Y) 

RETURN 

END 

(f) Subroutine ADJUST 



(g) Subroutine MULT 

$ I BFTC MULTD FULIST 

SUBROUTINE M UL T<M»N»K,A,B,C , IC , NA , MA , N& ,MB , NC , MC ) 

C THIS IS A GENERAL pURpOSE MATr t X MULTIPLICATION SUBROUTINE 

C 

DIMENSION A { NA , MA ) , B ( NB , MB ) , C { NC , MC } 

IF(IC) 10,20,30 
10 DO 15 1=1, N 
DO 15 J = 1 , K 
DO 15 L~ 1 > M 

15 C(I,J)=C(I , J ) -A ( I » L ) *B ( L , J ) 

RETURN 

20 DO 25 1 = 1, N 
DO 25 J= 1 , K 
C ( I , J ) = 0 • 

DO 25 L“ 1 » M 

25 C ( I » J ) -C (I , J ) +A ( I , L ) *B ( L , J ) 

RETURN 

30 DO 35 1=1, N 
DO 35 J= 1 , K 
DO 35 L= 1 , M 

35 C(I , J)=C( I, J)+A( I ,L)*B(L, J) 

RETURN 

END 
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$ I BMAP 

DOTSLV 

DOT , SDOT , DAD, 

ABD IL0G2 FUR USE WITH SOLVE(SLVIT) 

DOTS 

10 

❖ 


17 fEB. 1968 



DOTS 

20 

* 


R. J. HANSON 

, JPL. 


DOTS 

30 

❖ DOT AND FRIENDS ROUTINES FOR USE WITH SOLVE 

(SLVIT) 

DOTS 

40 


ENTRY 

DOT (n»A( 1) yMAf B( 1) , 

MB , C ) DOUBLE INNER PRODUCT 

DOTS 

50 


ENTRY 

SDOT (Ny A( 1 ) , MA » B( 1) 

, MB , C ) INNER PRODUCT 

DOTS 

60 


ENTRY 

IL0G2 (A) FLOATING 

POINT EXPONENT 


DOTS 

70 


ENTRY 

DA D (AyB) ADD WITH 

ROUND 


DOTS 

80 

❖ 





DOTS 

90 

SNAD 

MACRO 

M 

STD, COMPLEMENTING IF NECESSARY 

DOTS 

100 


CHS 




DOTS 

110 


ADD 

= 1B2 



DOTS 

120 


ALS 

18 



DOTS 

130 


STD 

M 



DOTS 

140 


ENDM 

SNAD 



DOTS 

150 

SAV 

MACRO 

A 



DOTS 

160 


SXA 

Ay 1 



DOTS 

170 


SXA 

A+1,2 



DOTS 

180 


SXA 

A+2,4 



DOTS 

190 


ENDM 

SAV,NOCRS 



DOTS 

200 

RETUR 

MACRO 

A 



DOTS 

210 


AXT 

**,1 



DOTS 

220 


AXT 

**,2 



DOTS 

230 


AXT 

❖❖ , 4 



DOTS 

240 


TRA 

1,4 



DOTS 

250 


ENDM 

RETUR, NOCRS 



DOTS 

260 

DOT 

SAV 

RET 



DOTS 

270 


STZ 

S 



DOTS 

280 


STZ 

S+l 



DOTS 

290 


CLA* 

8,4 



DOTS 

300 


LDO 

C + l 



DOTS 

310 


STO 

C 



DOTS 

320 


CL A* 

3,4 



DOTS 

330 


TZE 

NONE 

SKIP LOOP IF N - 

0 

DOTS 

340 


S TO 

N 



DOTS 

350 


CLA 

4,4 



DOTS 

360 


PAC 

, 1 

X 1 = - 1 BASE OF A) 


DOTS 

370 


CLA* 

5,4 

MA 


DOTS 

380 


SNAD 

MA 



DOTS 

390 


CLA 

6,4 



DOTS 

400 


PAC 

,2 

X2-- { BASE OF B) 


DOTS 

410 


CLA* 

7,4 



DOTS 

420 


SNAD 

MB 



DOTS 

430 


LXA 

N , 4 



DOTS 

44 0 

LOOP 

LDO 

0,1 

A ( I ) 


DOTS 

450 


FMP 

0,2 

B( I) 


DOTS 

46 0 


DEAD 

S 



DOTS 

470 


DST 

S 



DOTS 

480 

HA 

TXI 

* + 1 , 1,** 

(XI )={X1) +MA 


DOTS 

490 

MB 

TXI 

*+l , 2 , ** 

( X2 ) = ( X2 ) +MB 


DOTS 

500 


T I X 

LOOP, 4, 1 

END OF MAIN LOOP 


DOTS 

510 

NONE 

DEAD 

C 



DOTS 

520 


FRN 




DOTS 

530 

RET 

RETUR 




DOTS 

540 

* 





DOTS 

550 

SOOT 

SAV 

SRET 



DOTS 

560 


STZ 

S 



DOTS 

570 


CLA* 

8,4 



DOTS 

580 


STO 

C 



DOTS 

590 


CLA* 

3,4 



DOTS 

600 


TZE 

SNONE 



DOTS 

610 


STO 

N 



DOTS 

620 


CLA 

4,4 



DOTS 

630 


PAC 

,1 



DOTS 

640 


CLA* 

5,4 



DOTS 

650 


SNAD 

SMA 



DOTS 

660 


CLA 

6,4 



DOTS 

670 


PAC 

,2 



DOTS 

680 


CLA* 

7,4 



DOTS 

690 


snad 

SmB 



DOTS 

700 
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LXA 

N,4 

DOTS 710 

SLOOP 

LDQ 

°a 

D°TS 7 20 


FMP 

0,2 

DOTS 730 


FAD 

s 

DOTS 740 


S T0 

S 

DOTS 750 

SMA 

TXI 

*+1,1,** 

DOTS 760 

SMB 

TXI 

*+ 1 , 2 , ** 

DOTS 770 


T I X 

SLOOP, 4,1 

DOTS 780 

SNONE 

FAD 

C 

DOTS 790 


FRN 


DOTS 800 

SRET 

RETUR 


DOTS 810 

❖ 



DOTS 820 

ILOG2 

CAL* 

3,4 

DOTS 830 


ANA 

=0377000000000 

DOTS 840 


SUB 

=0200000000000 

DOTS 850 


ARS 

27 

DOTS 860 


TRA 

1,4 

DOTS 870 

* 



DOTS 880 

DAD 

CLA* 

3,4 

DOTS 890 


FAD* 

4,4 

DOTS 900 


FRN 


DOTS 910 


TRA 

1,4 

DOTS 920 

# 



DOTS 930 


EVEN 


DOTS 940 

C 

PZE 


DOTS 950 


PZE 


DOTS 960 

S 

PZE 


DOTS 970 


PZE 


DOTS 980 

N 

PZE 


DOTS 990 

T 

PZE 


DOT S 1000 


END 


D0TS102 0 


(h) Subroutine SOLVE 


$ I BFTC SLVEIT FUL I ST , DECK 



SLVE0010 



SUBROUTINE SOLVE IL , A, N, B , X , I N , EPS , I TMAX , IT,AA,S) 



SLVE0020 

CSOLVE 

LINEAR EQUATION SOLVER WITH ITERATIVE IMPROVEMENT 

VERSION II 

SLVE0030 

C 


SOLVES AX=B WHERE A IS NXN MATRIX AND B IS NX1 VECTUR 


SLVE0040 

C 






SLVE0050 

C 


IN = 




SLVE 0060 

c 



1 FOR FIRST ENTRY 



SLVEOO 7 O 

c 



2 FOR SUBSEQUENT ENTRIES WITH NEW B 



SLVE0080 

c 



3 TO RESTORE A AND B 



SLVE0090 

c 



4 IF FIRST ENTRY BUT AN INITIAL APPROXIMATION 

IS 

already 

SLVE0100 

c 



AVAILABLE IN X 



SLVE01 1 0 

c 



5 CONTINUE CALCULATING ITERATIVE IMPROVEMENT 

FOR 

THIS SYSTEM 

SLVE0120 

c 






SLVEOl 3 0 

c 


EPS 

AND I TMAX ARE PARAMETERS IN THE ITERATION 



SLVE0140 

c 


I T = 




SLVE015 0 

c 



-1 IF A IS SINGULAR 



SLVE0160 

c 



0 IF NOT CONVERGENT 



SLVEOl 70 

c 



NUMBER OF ITERATIONS IF CONVERGENT 



SLVE0180 

c 


CALLS MAP SUBROUTINES IL0G2, DOT, SDOT AND DAD 



SLVE0190 

c 






SLVE0200 



DIMENSION A(L,L) ,B(L) ,X(L) ,AA(N,N) ,S(1) 



SLVE0210 



MA = 

L 



SLVE0220 

c 


MA 

MUST = DECLARED DIMENSION OF SYSTEM 



SLVE023 0 



Jl = 

N 



SLVE0240 



J2 = 

Jl+N 



SLVE0250 



J 3= 

J2+N 



SLVE0260 



GO 

TO ( 10,230,390,10,290) ,IN 



SLVE0270 


10 

NM1 

=N — 1 



SLVE0280 



NP1 

= N+ 1 



SLVE0290 

c 






SLVE0300 

c 


EQUILIBRATION 



SLVE031 0 

c 






SLVE0320 



DO 

40 I = 1 , N 



SLVE0330 




KT0P = IL0G2( At I , 1 ) ) 



SLVE0340 




DO 20 J = 2 » N 



SLVE035 0 


20 

KT0P=MAX0{KT0P,IL0G2<A(I,J) } ) 



SLVE0360 



N2= 

J2+I 



SLVE0370 



S(N2) = 2 . 0 **(-KT 0 P ) 



SLVE0380 




DO 30 J= 1 » N 



SLVE0390 


60 
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30 

All , J)=A< I , J)*S(N2) 



SLVE0400 


40 

CONTINUE 



SlvE°410 

c 





SLVE0420 

c 


SAVE EQUILIBRATED DATA 



SLVE0430 

c 





SLVE0440 



DO 50 1=1, N 



SLVE0450 



DO 50 J=1,N 



SLVE0460 


50 

AA ( I » J ) = A ( I , J) 



SLVE0470 

c 





SLVE0480 

c 


GAUSSIAN ELIMINATION WITH PARTIAL 

PIVOTING 


SLVE0490 

c 





SLVE0500 



DO 170 M=1,NM1 



SLVE0510 



N3=J3+M 



SLVE0520 



TOP=ABS ( A ( M , M ) ) 



SLVE0530 



I M AX=M 



SLVE0540 



DO 70 I =M , N 



SLVE0550 



I F { TOP-ABS ( A ( I ,M) }) 60,70,70 



SLVE0560 


60 

TOP=ABSI A( I ,M) ) 



SLVE0570 



I MAX= I 



SLVE0580 


70 

continue 



SLVE0590 



IF ( TOP ) 90 , 80 ,90 



SLVE0600 


80 

I T=-l 



SLVE0610 

c 


^SINGULAR* 



SLVE0620 



return 



SLVE0630 


90 

S { N3 ) = I MAX 



SLVE0640 


100 

I F ( IMA X-M )1 30, 130, 110 



SLVE0650 


HO 

DO 120 J= 1 ,N 



SLVE0660 



T EMP=A ( M , J ) 



SLVE0670 



AIM, J) = A( I MAX , J ) 



SLVE0680 


120 

A ( I MAX , J ) =T EMP 



SLVE0690 


130 

MP 1 = M+1 



SLVE0700 



DO 160 I =MP 1 , N 



SLVE0710 



EM= A ( I , M ) / A ( M, M ) 



SLVE0720 



A ( I , M ) =EM 



SLVE0730 



I F ( EM ) 1 40 , 1 60 , 1 40 



SLVE0740 


140 

DO 150 J = MP 1 » N 



SLVE0750 



I F ( I L0G2 ( A ( M , J ) ) + ILoG2(EM) .LT. 

-54) GO TO 

150 

SLVE0760 



A ( I , J ) = A ( I » J)-A<M, J)*EM 



SLVE0770 

c 

150 

A( I , J)=A( I , J)-A(M, J)*EM 



SLVE0780 


150 

CONTINUE 



SLVE0790 


160 

continue 



SLVE0800 


170 

CONTINUE 



SLVE0810 



N4=N*4 



SLVE0820 



S (,N4 ) =N 



SLVE0830 



IF { A ( N , N ) ) 190 , 180 ,190 



SLVE0840 


180 

I T=“l 



SLVE0850 



RETURN 



SLVE0860 


190 

CONTINUE 



SLVE0870 

c 


STORAGE FOR A NOW CONTAINS TRIANGULAR L AND U 

SO THAT ( L+I)*U=A 

SLVE0880 

c 





SLVE0890 

c 


DUPLICATE INTERCHANGES IN DATA 



SLVE0900 

c 





SLVE091 0 



DO 220 1=1, N 



SLVE0920 



N3= J3+ 1 



SLVE0930 



I P=S ( N3 ) 



SLVE0940 



IF ( I - IP ) 200 , 220 i 200 



SLVE0950 


200 

DO 210 J= 1 , N 



SL V E 0960 



TEMP = AAU , J) 



SLVE0970 



AA ( I » J ) = AA ( IP, J) 



SL VE 09 80 


210 

AA { I P , J ) = T EMP 



SL VE 099 0 


220 

CONTINUE 



SLVE1000 

c 





SLVE1010 

c 


PROCESS RIGHT HAND SIDE 



SLVE1020 

c 





SL VE 1 03 0 


230 

DO 240 1=1, N 



SLVE1040 



N2= J2+ 1 



SL VE 1 05 0 


240 

B ( I )«B C I ) *S (N2 ) 



SLVE1060 



DO 250 1=1, NM1 



SL VE 107 0 



N3=J3+I 



SLVE1080 



I P= S { N 3 ) 



SL VE 109 0 
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Printout D-2 (contd) 




T EMP = B ( I ) 

SLVE1100 



B ( I >=B<IP) 

slveUIO 



B ( Ip)=TEMp 

SLVEU20 


250 

CONTINUE 

SLVE1130 



IF ( IN-4) 260,290,290 

SLVE1140 

c 



SLVE1150 

c 


BYPASS INITIAL APPROXIMATION CALCULATION 

SLVE1 160 

c 


IN GREATER THAN 4 CANNOT OCCUR AT THIS POINT 

SLVE1170 

c 



SLVE11 80 

c 



SLVE1190 

c 


SOLVE FOR FIRST APPROXIMATION TO X 

S L VE 1200 

c 



SLVE1210 


260 

DO 270 I = 1 » N 

SLVE1220 



Nl= J 1 +1 

SLVE1230 


270 

S<N1) = -SD0T( 1-1, A< I ,1) ,MA,S( Jl + 1) , 1,-B( I } ) 

SLVE1240 



DO 280 K= 1 , N 

SLVE1250 



I =N P 1-K 

SLVE1260 



Nl-Jl+I 

SLVE1270 


280 

X ( I )=-SDOT ( N — I » A( I , 1+1) ,MA,X( 1 + 1) ,lt-S(Nl) )/A( I ,1 ) 

SL VE 1280 

c 



SLVE1290 

c 


iterative improvement 

SLVE13 00 

c 



SLVE1310 


290 

1 F { ITMAX) 370,370,300 

SLVE1320 


300 

TGP=0 .0 

SLVE1330 



DO 310 I = 1 , N 

SL VE 134 0 


310 

TOP = AMAX1(TOP,ABS{X{ I ) ) ) 

SLVE1350 



EPSX=EP$*TOP 

SLVE1360 



DO 360 IT=1, ITMAX 

SLVE1370 

c 


FIND RESIDUALS 

SLVE1380 



DO 320 1=1, N 

SLVE1390 


320 

S ( I ) =— DOT ( N » A A ( I i 1 ) * N,X( 1) , 1,-B( I ) ) 

SL VE 1400 

c 


FIND INCREMENT 

SLVE1410 



DO 330 1=1, N 

SLVE1420 



N 1= Jl+I 

SLVE1430 


330 

S(N1)=-$D0T ( I -1 , A ( I ,1) , MA, S( Jl+1 ) , 1,-S(I ) ) 

$L VE 1440 



DO 340 K= 1 ,N 

SLVE1450 



I =NP1-K 

SLVE1460 



Nl= Jl+I 

SLVE1470 


340 

S ( I } = -SDOT(N-I ,A( I ,1 + 1) ,MA,S( 1 + 1 ) t lt-SINl) )/A{I ,1 ) 

SLVE1480 

c 


INCREMENT AND TEST CONVERGENCE 

SL VE 149 0 



1 UP=0*0 

SLVE1500 



DO 350 1=1, N 

SLVE1510 



TEMP=X(I ) 

SLVE1520 



X(I )=DAD(X(I),SU)} 

SLVE1530 



DELX = AB$ ( X 1 1 ) -TEMP ) 

SLVE1540 



TOP = AMAXUTOP,DELX) 

SLVE1550 


350 

continue 

SLVE1560 



IF (TOP-EPSX) 380,380,360 

SLVE157 0 


360 

CONTINUE 

SLV'E 15 80 


370 

I T = 0 

SLVE1590 


380 

RETURN 

SLVE1600 

c 



SLVE161 0 

c 


RESTORE A AND B 

SLVE1620 

c 



SLVE1630 


390 

CONTINUE 

SLVE1640 



DO 420 K=1,N 

SLVE1650 



I =N P 1-K 

SLVE1660 



N3= J3+ I 

SL VE 16 70 



IP=S(N3) 

SLVE1680 



IF U-IP >400,420,400 

SLVE1690 


400 

T EM P=B ( I ) 

SLVE1700 



B ( I ) = B( IP) 

SLVE1710 



B(IP)=TEMP 

SLVE 1720 



DO 410 J=1,N 

SLVE1730 



TEMP=AA ( I , J) 

SLVE 1740 



AAII , J)=AA(IP, J) 

SLVE 1 750 


410 

AA ( I P , J ) =TEMP 

SLVE 1760 


420 

CONTINUE 

SLVE1770 



DO 430 1=1, N 

SLVE 1780 



N2= J2+ 1 

SLVE 1 790 
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o o o o o n 


Printout D-2 (contd) 


B(I ) = 8 ( I )/S(N2) 

DU 430 J=1,N 
A(I,J)=AAU,J)/$(N2) 

430 CONTINUE 
RETURN 
END 

$1 BFTC SLVIN. DECK 

SUBROUTINE SLVINV (NDIM,N,A,B,IT,S) 

THE DIMENSION OF S MUST BE AT LEAST N**2+5*N 


SLV 

DIMENSION B(NDIM,NDIM) ,S(N) SLV 

ITMAX » 10 SLV 

IN = 1 SLV 

JJ = 4*N SLV 

DO 20 J = 1 » N SLV 

DO 10 I ~ 1,N 
K - JJ + I 
10 S(K) = 0.0 
K = JJ + J 
S ( K ) = 1.0 


CALL SOLVE ( NDI M , A , N , S ( 4*N+ 1 ) ,B( 1, J ) » IN , 7 . 0E-9 > IT MAX, IT, Si 
IF (lT.EQ.t-l) ) GO TO 30 
IF (IT.EQ.O) CALL SOLVE ( ND I M , A , Nt S ( 4*N+ 1 ) , B ( 1 , J ) , 5 , 7 .0 E-9 , IT MAX , 
1 IT,S(5*N+1 ) ,S) 

20 IN = 2 
CONTINUE 


RESTORE THE MATRIX A 

30 CALL S0LVE(NDIM,A,NfS<4*N+l),B(ltJ), 
RETURN 
ENO 


SLVE1800 
SLVE 1810 
SLVE1820 
SLVE1830 
SLVE 1 840 
SLVE 1850 
SLV I 001 0 
SLV 1 0020 
SLV I 003 0 
SLV I 0040 
10050 
10060 
1007 0 
10080 
10090 
10100 
SLV I 01 1 0 
SLV 1 0 1 20 
SLV I 0130 
SLV 1 0140 
SLV I 0 1 5 0 
5*N+1) » S ) SL V 1 0160 
SLVI0170 
SLV 10180 
SLV I 0 1 90 
SLV 1 0200 
SLV I 02 1 0 
SLV I 0220 
SLV 1 0230 
SLV 1 0240 

3,7.0E-9, ITMAX, IT,S(5*N+1) ,S)SLVI0250 

SLV 1 0260 
SLV I 0270 
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Nomenclature 


Coordinates and Related Quantities 


x,y,z 

Cartesian coordinates (basic coordinate 
system with the origin located at the 
source of the incident electromagnetic 
field) 

P 

field point (the point at which the scat- 
tered or resultant field is evaluated) 

R 

the vector from the origin to P 

R, 0,5- 

spherical coordinates of P 

S 

the scattering surface 

P 

the vector from the origin to a point 
on S 

p, $, <j> 

spherical coordinates of a point on S 

p(6, </») 

a function describing S 

AS mn 

an incremental area on S 

r 

the distance from P to a point on S 

n 

a unit vector normal to S 

i 

a unit vector in the direction indicated 
by a coordinate subscript 

P 

feed offset angle 

V 

a region of space 

Electromagnetic Quantities 

E 

electric field 

H 

magnetic field 

Ei,H; 

incident field 

E s , H s 

scattered field 

Ey 

total resultant field 

Hp, H s , H<j, 

components of H* with radial variation 
factored out 

Ee,E$ 

components of Ej with radial variation 
factored out 

Eq, e® 

components of E s with radial variation 
factored out 

K 

surface currents 

AK 

difference between true currents and 
physical optics currents 

k = a>(£//.) Vs — 2nr/X 

the propagation constant 


6 electric permittivity 
\ wavelength 
fx magnetic permeability 
co angular frequency 


Functions and Constants Used in Derivations 


I 


Ie, 1$ 


Al 


mn 


ajfi), hn{0) 

& e 0 mn > h * mn 


a vector function related to E s 
components of I 

incremental contribution to I from AS wm 
Fourier components of surface data 

coefficients related to the spherical- 
wave coefficients 


0 n l spherical- wave coefficients 

b n b^ mn ^ 

n ?l =n« m J 

> spherical-wave functions 
y the path-length function 


F a vector function related to H t and p 


a m , n? b wm > c mn interpolation coefficients for the func- 
tion F 


A»SW=A((9)\ 

B m ° 0 (d)=B(e)) 


interpolation coefficients for the func- 
tion y 

mth even or odd Fourier component of 
incident electric field 


Miscellaneous Symbols 


j (-ir 

V Laplacian operator 
fi n (kp) spherical Bessel function 
P”* (cos 0) associated Legendre function 
f focal length of a paraboloid 


D physical diameter of source (also used 
as diameter of a scattering surface) 

Fj (0), G™ (0) functions related to the associated 
Legendre function 
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